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ABSTRACT 


The  continuous  integrals  and  integral  equations  which 
form  the  theory  of  Nearfield  Acoustic  Holography  for  planar 
and  odd-shaped  source  boundary  surfaces  are  reviewed,  and 
the  approximations  and  assumptions  necessary  to  reduce  these 
equations  to  a  set  of  finite  and  discrete  operations 
suitable  for  computation  are  developed.  These  equations 
represent  the  solution  of  the  Helmholtz  equation  with 
specified  boundary  conditions  by  Green's  function  methods. 

In  analyzing  the  reduction,  to  discrete  form,  of  the 
propagation  of  planar  holograms  (a  record  of  the  radiated 
field  over  a  plane  above  the  boundary  plane)  ,  four  new 
methods  for  representing  the  Green's  functions  numerically 
are  developed.  The  relationship  of  two  earlier  Green's 
function  representation  methods  to  these  new  forms  is 
established.  The  results  of  numerical  testing  of  these  six 
representations  are  presented.  These  results  demonstrate 
the  effectiveness  of  the  new  representations. 

Two  computational  methods  for  reconstructing  planar 
source  boundary  fields  from  planar  holograms  are  developed. 
The  first  method  is  an  approximation  of  the  continuous 
solution  method  which  the  convolution  theorem  of  Fourier 


t 


p.  i 

Transforms  provides.  In  its  simplest  form,  this  method  has 
a  high  sensitivity  to  noise;  it  is  shown  how  this  problem  is 
partially  alleviated  by  the  use  of  a  high  spatial  frequency 
filter.  Numerical  test  results  indicate  the  effectiveness 
of  this  method  for  very  short  distances  between  hologram  and 
source  boundary,  a  condition  often  satisfied  in  experiment. 
The  new  methods  for  representing  the  Green’s  function, 
developed  for  propagating  planar  holograms,  are  shown  to  be 
inappropriate  with  this  reconstruction  technique. 

As  an  alternative  to  this  reconstuction  method,  a 
conjugate  gradient  descent  method  is  developed  based  on  the 
finite  and  discrete  propagation  method  discussed v  This 
technique  does  benefit  from  the  new  Green's  function 
representations.  Although,  this  reconstruction  method 
requires  more  time  than  the  first  method  presented,  it  is 
relatively  insensitve  to  noise,  and  it  extends  the  feasible 
reconstruction  range,  or  distance  between  hologram  and 
boundary . 

'.The  reduction  to  finite  and  discrete  form,  by  a  Finite 
Element  technique,  of  the  relationship  between  the  Dirichlet 
and  Neumann  boundary  conditions  for  an  odd  shaped  surface  is 
reviewed.  To  develop  a  unique  relationship,  assuming 
radiation  conditions  apply,  a  knowledge  of  the  boundary 
surface's  characteristic  frequencies  is  essential,  and  a 
method  for  detecting  these  frequencies  for  an  odd-shaped 


boundary  is  presented.  A  comparison  of  the  characteristic 
frequencies  predicted  by  this  method  and  those  given  by 
theory  is  presented  for  a  spherical  surface. 


A  technique  for  reconstructing  the  odd-shaped  surface 
boundary  conditions  from  a  hologram  of  general  two-dimen¬ 
sional  shape  is  proposed.  Results  from  a  numerical  study 

> 

which  support  this  technique  are  presented.  Two  cases  are 
considered  in  this  study:  reconstruction  of  a  spherical 
boundary  from  a  planar  hologram,  and  reconstruction  of  a 
spherical  boundary  from  a  concentric,  spherical  boundary. 
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Chapter  I. 


INTRODUCTION  TO  NEARFIELD  ACOUSTIC  HOLOGRAPHY  (NAH) 

The  works  presented  in  this  thesis  are  refinements  and 
extensions  to  the  comptational  theory  and  practice  of  an 
exciting  new  imaging  technique,  conceived  and  developed  at 

Penn  State,  called  Nearfield  Acoustic  Holography.1  The 
computational  techniques  of  Nearfield  Acoustic  Holography 
are  of  particular  interest  to  Physists,  as  they  represent 
the  development  and  application  of  classical  field  theory  to 
large  scale  and  geometrically  complex  problems.  The 
solution  of  such  problems  has  become  feasible  only  recently 
as  a  result  of  the  explosive  development  in  computer 
facilities  in  the  last  two  decades.  Quite  naturally,  it  is 
this  author's  hope  that  the  methods  explored  and  developed 
here  will  continue  to  find  application  in  other  areas  of 
Physics,  but  for  now  this  exposition  must  be  limited  to  the 
specific  problems  encountered  in  Nearfield  Acoustic 
Holography . 

In  the  first  section  of  this  chapter,  the  basic 
concepts  of  holography,  acoustic  holography,  and  Nearfield 
Acoustic  Holography  are  presented  in  turn.  The  major 
impetus  for  developing  NAH  was  its  potential  as  an 

experimental  technique,  and  defining  the  experimental 


2)  the  subsequent  display  of  any  region  of  the  original 
three-dimensional  field  based  on  the  single 
two-dimensional  record  of  the  field  obtained  in  part  1 . 
The  record  of  this  surface  field  information  is  refered  to 
as  the  hologram. 

The  ideal  surface,  over  which  the  information  is 

recorded,  the  hologram  surface,  divides  space  into  a  region 

with  sources  and  a  region  without  sources.  The  ideal 

hologram  contains  all  information  about  the  field  over  this 

entire  surface.  Of  course,  in  practice,  these  ideals  are 

neither  possible  nor  necessary,  and,  hence,  the  requirement 

to  specify  the  field  uniquely  is  relaxed  to  a  practical 

level.  It  is  the  particular  application  which  imposes  the 

minimum  information  requirements.  For  example,  when  viewing 
* 

an  optical  hologram  of  a  set  of  chess  pieces,  one  would  not 
expect  to  resolve  micron  size  chips  in  the  surface  of  a 
pawn,  even  though  such  information  is,  in  theory,  carried  by 
the  original  field.  More  specific  requirements  as  they, 
apply  to  NAH,  are  discussed  in  later  chapters. 

A. 2.  Acoustic  Holography 


of  particle  velocity,  and  the  three  components  of  acoustic 
intensity  can  be  determined  at  any  point  in  the  source-free 
region  of  the  acoustic  medium.  NAH  essentially  meets  this 
goal,  although  acoustic  holography  techniques  of  the  past 
fell  quite  short  of  this  goal  in  experimental  theory  and 
practice . 

Previous  forms  of  acoustic  holography  were  largely 
based  on  the  experimental  methods  of  optical  holography,  and 
many  required  the  use  of  lasers  to  redisplay  the  acoustic 
fi  — 7 

field.  These  methods  recorded  either  a  component  of 

particle  velocity  or  the  acoustic  pressure  over  the  hologram 
surface.  Since  these  methods  depended  on  optical  inter¬ 
ference  phenomena  to  redisplay  the  holographically  recorded 
field  component,  subsequent  observation  of  the  reproduced 
field  had  to  be  accomplished  by  a  phase-  insensitive  optical 
device.  The  resulting  lack  of  phase  information  precluded 
the  accurate  determination  of  the  remaining  acoustic  field 
variables,  including  the  acoustic  intensity  vector  field. 

Another  problem  with  previous  acoustic  holography 
techniques  was  resolution.  This  problem  resulted  from 
developing  acoustic  holography  from  experimental  optical 
holography  whose  theory  and  success  is,  for  the  most  part, 
based  on  the  assumed  existence  of  a  purely  propagating  wave 
field.  To  insure  the  validity  of  this  assumption,  a 
hologram  must  be  located  in  the  farfield  of  any  source. 


Since  they  were  based  on  experimental  optical  methods, 
earlier  forms  of  acoustic  holography  were  similarly 
restricted.  This  restriction  is  of  little  consequence  in 
optical  holography,  although  it  imposes  a  wavelength 
resolution  limit  on  the  restored  fields.  However,  if  this 
resolution  limit  was  physically  inherent  in  all  holography 
theory,  acoustic  holography  would  be  useless  for  many 
problems  in  the  audio  range,  where  wavelengths  are  typically 
larger  than  their  source.  Fortunately  this  limitation  is 
only  a  practical  limitation  occurring  in  most  forms  of 

optical  holography,8  and  it  may  be  viewed  as  artificial  in 
acoustic  holography. 


A. 3.  Nearfield  Acoustic  Holography 

Nearfield  Acoustic  Holography  is  an  acoustic 
holography  technique  which  comes  much  closer  to  the  ideal 
holography  described  above.  As  its  name  declares,  this 
method  is  carried  out  in  the  nearfield  of  radiating  sources, 
and  much  of  its  merit  comes  from  the  ability  of  this 
technique  to  accurately  process  the  non-propagating  wave 
field  information  present  in  the  nearfield.  The  resolution 
of  this  technique  is  not  limited  by  the  wavelength. 


To  understand  how  this  increase  in  resolution  is 
obtained,  and  why  former  techniques  could  not  achieve  this 
resolution,  it  is  expedient  to  separate  the  possible  finite 
linear  acoustic  wave  solutions,  in  a  semi-infinite  space, 
into  two  catagories:  propagating  waves  which  transport 
energy  out  to  infinity,  and  evanescent  waves  which  are 

O 

characterized  by  their  rapid  decay.  As  discussed  in 

chapter  II  and  the  original  NAH  paper1  for  planar  bound¬ 
aries,  this  separation  is  performed  explicitly  by  the 
Fourier  transform,  and  the  resulting  evanescent  waves 
manifestly  decay  exponentially.  Information  about  details 
of  a  source  smaller  than  the  characteristic  acoustic  wave¬ 
length  can  only  be  obtained  by  measuring  the  strength  of 
these  evanescent  waves,  restoring  them  to  their  original 
strength  at  the  source,  and  finally  superimposing  these 

waves  with  the  readily  restored  propagating  waves.1  It  is 
not  possible  to  restore  these  decayed  waves,  a  process  which 
requires  selective  enhancement,  in  the  older  forms  of 
acoustic  holography  which  rely  on  optical  wave  propagation 
to  achieve  a  representative  optical  restoration  of  the 
acoustic  field. 

NAH  does  not  require  such  a  restoration  technique  due 
to  its  unique  process  for  forming  the  hologram;  the  hologram 
is  an  electronically  recorded  sample,  in  digital  form,  of 
the  field  at  a  large  number  of  points  over  the  hologram 
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surface.  With  this  type  of  hologram,  numerical  restoration 
of  evanescent  waves  is  relatively  easy,  and  phase 
information  is  also  readily  obtained  and  manipulated 
numerically .  The  numerically  restored  fields  are  displayed 
with  computer  graphics. 

At  this  point,  an  examination  of  the  experimental 
hardware  requirements  of  NAH  will  define  the  particulars  of 
the  overall  technique. 

B.  Experimental  requirements 

The  experimental  implementation  of  NAH,  such  as  the 

prototype  system  at  Penn  State,9  consists  of  the 
successful  integration  of  an  extensive  system  of  electronic 
equipment.  An  appreciation  of  the  capabilities  and  limita¬ 
tions  of  such  a  system,  gained  by  examining  its  experimental 
requirements,  will  help  clarify  the  restrictions  and 
approximations  introduced  in  the  remaining  chapters,  where 
computational  NAH  methods  are  developed  largely  with 
experimental  systems  in  mind. 

The  experimental  equipment  may  be  divided  into  two 
subsystems:  the  transducers  and  the  computer  facilities. 

The  requirements  and  tasks  for  each  subsystem  are  outlined 


below . 


B.l.  Transducers 


Inorder  to  avoid  aliasing,  the  transducer  subsystem  in 
a  typical  NAH  set-up  must  supply  the  computer  subsystem,  at 

a  rate  at  least  twice  the  maximum  acoustic  frequency,10  with 
a  digital  representation  of  the  acoustic  pressure  at  a  each 
of  the  points  which  define  the  hologram  surface.  This  sub¬ 
system  must  also  supply  the  exact  location  of  the  measure¬ 
ment  points.  Two  methods  have  been  employed  to  locate  the 
hologram  measurement  points,  and  these  two  methods  also  have 
different  equipment  requirements  and  data  conversion  limita¬ 
tions  . 

A  16x16  square  array  of  microphones  forms  the  heart  of 
the  prototype  transducer  system  at  Penn  State.  Each 
microphone  is  fixed  in  position  in  the  array,  and  the  entire 
array  is  positioned  by  synchronous  motors  under  computer 

control.11  With  this  system  it  is  easy  to  determine  the 
relative  and  absolute  positions  of  each  microphone.  Each 
microphone  is  equipped  with  a  preamplifier,  noise  filter, 

and  an  analog  switch.12  As  part  of  a  relatively  complex 
high  speed  switching  system,  the  analog  switch  selectively 
connects  the  microphone  to  a  shared  signal  line,  and  this 
signal  is  then  routed  to  a  high  speed  analog  to  digital 
1  3 

converter.  The  advantages  of  such  a  system  are  the 


speed  with  which  single  frequency  sources  can  be  analyzed 
and,  due  to  the  nearly  simultaneous  sampling  of  the 
microphones  achieved  with  the  switching  system,  the 
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possibility  of  analyzing  relatively  wideband  sources. 

Some  of  the  complexity  and  speed  of  a  microphone  array 
system  can  be  traded  off  in  the  study  of  single  frequency 
sources,  or  highly  reproducible  wideband  sources,  by  using 
an  articulated  boom  mounted  microphone.  Besides  the  loss  in 
speed,  the  use  of  this  type  of  system  requires  a  more 
complex  system  for  specifying  the  microphone  location.  On 
the  positive  side,  with  a  single  microphone,  phase  sensitive 
detection  methods  can  be  employed  to  reduce  the  single  to 

noise  ratio  in  single  frequency  studies.14  Also,  an 
articulated  boom  mounted  microphone  allows  the  specification 
hologram  surfaces  which  are  not  flat,  and  this  may  be  of 
great  advantage  in  the  study  of  odd-shaped  sources  as 
discussed  in  chapters  V  and  VI .  Boom  mounted  systems  have 
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been  employed  by  a  number  of  investigators.  '  In  this 
type  of  system,  the  final  output  is  either  a  near  real-time 
digital  representation  of  the  microphone  signal,  which  is 
collected  and  processed  by  the  computer  system  or,  for 
single  frequency  work  with  phase-sensitve  detection,  a 
digital  representation  of  the  phase  ( e . g ., relative  to  the 
source  drive  signal)  and  amplitude  of  the  microphone  signal. 


The  computer  system  in  a  NAH  experiment  serves  to 
store,  process,  and  display  all  collected  and  generated 
data,  and  it  is  easiest  to  review  the  system  requirements  by 
describing  the  Penn  State  system  and  the  role  its  components 
play.  The  Penn  State  system  is  composed  of  a  minicomputer 
with  the  following  attached  components:  a  peripheral  port 
connection  with  the  analog  to  digital  system,  hard  disk  and 
floppy  disk  storage  devices,  an  array  processor,  a  vector 
graphics  display  monitor,  and  a  pen  plotter.  The  role  of 
each  of  these  devices  is  seen  by  following  the  typical 
sequence  of  steps  in  a  NAH  experiment. 

A  digital  sample  of  the  signal  at  each  microphone 
location,  as  a  function  of  time,  is  collected  by  the 
computer  from  the  analog  to  digital  conversion  system.  Each 
of  these  time  sequences  is  passed  to  the  array  processor 
where  they  are  transformed  into  a  frequency  domain 
representation.  An  important  step  in  this  part  of  the 
process  is  calibration.  In  a  multiple  microphone  system, 
calibration  information  is  stored  for  each  microphone 
providing  its  relative  phase  and  amplitude  response.  The 
measured  response  of  each  microphone  is  normalized  with  this 
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calibration  data.  Taken  at  a  single  frequency,  the  set 
of  complex  numbers,  representing  the  acoustic  pressure 


amplitude  and  phase  at  each  microphone,  defines  an  acoustic 
hologram.  Any  number  of  these  single  frequency  acoustic 
holograms  are  stored  on  magnetic  disks.  From  one  of  these 
holograms,  a  digital  representation  of  the  original  acoustic 
field  at  any  number  of  spatial  locations  is  generated  by 
numerical  methods  (these  methods  being  the  point  of  this 
thesis).  Since  an  enormous  amount  of  data  is  normally 
generated,  it  is  usually  displayed  on  one  of  the  graphic 
output  devices.  This  graphic  output  typically  takes  the 
form  of  a  hidden-line  plot  of  an  acoustic  variable  versus 
two  spatial  dimensions  as  in  the  example  of  figure  1.1,  or 
it  may  take  the  form  of  a  vector  intensity  map  as  in  the 
example  of  figure  1.2.  The  time  dependence  of  the  acoustic 
parameters  be  displayed  with  motion  pictures  on  the  display 
1  8 

monitor.  This  type  of  display  is  particularly  instructive 
in  the  study  of  source  boundaries,  as  the  source  structure's 
complicated  dynamic  response  may  be  veiwed  in  slow-motion. 

C.  Applications 


The  direct  applications  of  NAH  are  numerous  due  to  the 
wealth  of  information  about  an  acoustic  field  it  can  provide 
from  a  relatively  small  amount  of  input  data.  There  is  much 
interest  in  acoustic  intensity  measurement  today,  and  the 
holographic  method  can  quickly  provide  a  detailed  mapping  of 
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this  energy  flow  as  can  no  other  technique,  such  as  the  "two 

1  Q 

microphone  technique",  provide  in  a  reasonable  time.  An 

20 

advantage  which  NAH  has  over  dual  or  multiple  microphone 
texhniques  is  that,  the  intensity  determined  at  a  specific 
point  by  the  holographic  process  is  dependent  on  a  large 
number  of  microphone  measurements  (typically  4096  with  the 
prototype  system  at  Penn  State),  whereas  the  intensity 
determined  with  other  techniques  is  dependent  only  a  few 
measurements.  This  integration  of  a  large  number  of 
measurements  can  be  expected  to  improve  the  signal  to  noise 
ratio,  relative  to  other  field  mapping  techniques. 

Its  speed  and  accuracy  make  NAH  ideally  suited  for  the  study 
of  sound  production  equipment,  and  prototype  systems  for 

studying  musical  instruments  and  loudspeakers  are  already  in 
use.15'21 

NAH  also  provides  a  fast  non-contact  method  for 
determining  the  normal  component  of  surface  velocity  of  a 
vibrating  source.  A  particularly  interesting  application 
of  this  capability  is  in  the  proposed  experimental 
observation  of  an  acoustic  analogue  of  2-D  Anderson 
2  2 

Localization,  where  the  technique  will  be  used  to  observe 


the  vibrational  modes  of  a  randomly  nonuniform  elastic 

o  o 

sheet.  From  the  holographically  determined  normal 
component  of  velocity,  the  flow  of  energy  within  a  plate  may 

2  4 

also  be  studied. 

Any  area  where  quick  and  accurate  determination  of  the 
acoustic  field  above  a  plane  or  other  level  surface  in  a 
separable  coordinate  system  will  benefit  from  the 

Q 

application  of  NAH  in  its  present  state  of  development. 

With  the  experimental  implementation  of  the  extensions  of 
NAH  to  odd-shaped  sources,  and  the  refinement  of  the 
computational  methods  of  NAH  presented  in  this  thesis,  the 
list  of  practical  applications  should  grow. 


Chapter  II. 


FUNDAMENTAL  EQUATIONS  OF  NEARFIELD  ACOUSTIC  HOLOGRAPHY 

The  numerical  methods  of  Nearfield  Acoustic  Holography 
presented  in  the  following  chapters  are  approximations  to  a 
set  of  continuous  integrals  and  integral  equations  relating 
the  acoustic  field  to  appropriate  boundary  conditions.  It  is 
the  purpose  of  this  chapter  to  introduce  these  equations  and 
to  clarify  the  objectives  of  this  work  in  terms  of  these 
equations . 

The  first  section  of  this  chapter  very  briefly  reviews: 
the  formulation  of  Newton's  Law  for  continuous  media, 

Euler's  equation,  which  provides  a  direct  link  between  the 
acoustic  pressure  and  particle  velocity  fields;  and  the 
linear  approximation  of  the  general  equations  of  non-viscous 
hydrodynamics,  which  yields  the  Wave  Equation  and,  for 
single  frequencies,  the  Helmholtz  equation. 

The  continuous  solution  of  the  Helmholtz  Equation  with 
boundary  conditions  is  accomplished  by  a  Green's  function 
technique.  The  basic  ideas  behind  these  techniques  are 
presented  in  section  B.  The  Helmholtz  Integrals  and 
Rayleigh's  Integrals  are  reviewed  in  this  section,  as  they 
follow  from  the  solution  of  the  Helmholtz  Equation  by 
Green's  function  methods. 


The  acoustic  field  variables  of  pressure  and  three 
components  of  particle  velocity  are  of  -l  ief  interest  in 


Nearfield  Acoustic  Holography,  and  techniques  are  developed 
for  determining  these  fields  as  generated  from  single 
frequency  sources.  Note  that  the  term  "source"  refers  to  a 
boundary  condition  in  this  thesis,  since  very  few  problems 
in  Acoustics  involve  any  true  volume  sources,  and  no  volume 
sources  are  considered  here.  Also,  considering  only  single 
frequency  sources  is  not  restrictive  as  the  extension  of 
these  techniques  to  multiple  frequencies  is  readily 
accomplished  by  superposition  of  a  number  of  single 
frequency  solution  fields. 

A . 1  Euler's  Equation 


Newton's  Second  Law  in  a  form  suitable  for  use  in 
problems  of  continuous  media  provides  the  direct  link 
between  the  Acoustic  pressure  field  and  the  three  components 
of  the  particle  velocity  field.  By  applying  Newton's  Second 
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Law  to  an  infinitesimal  volume  of  inviscid  fluid  with 
density  p  acted  upon  only  by  the  pressure  field  p  in  which 
it  is  immersed,  it  is  straightforward  to  show  that  the 
acceleration  of  this  volume  is  given  as 

t( t.t)  =  =  -~vp(r , t )  (2.1). 


This  is  Euler's  Equation.25 

All  solutions  to  problems  in  this  work  are  formulated 
in  terms  of  single  temporal  frequency  components.  A  general 
explicit  function  of  time,  f(t),  is  transformed  to  one 
explicitly  of  frequency,  F(u),  by  the  Fourier  transform 
defined  by 


F  ( w) 


’ 

f (t)eiutdt 

-oo 


(2.2)  , 


and  its  inverse 


f(t)  -  fa 


+  oo 


F(«)e  iutdo 


(2.3)  . 


By  expressing  the  velocity  and  pressure  in  equation  (2.1)  in 

terms  of  their  frequency  domain  representations,  v(r,u)  and 

p(r,u),  with  integrals  of  the  form  of  equation  (2.3), 

Euler's  Equation,  equation  (2.1),  can  be  written  as 

v(r,w)  =  ~vp(r,u)  (2.4). 


This  equation  allows  one  to  calculate  the  particle  velocity 
field  once  the  acustic  pressure  field  p(r,u)  has  been  found. 
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A. 2  The  Wave  and  the  Helmholtz  Equations 


The  solution  of  the  general  equations  of  inviscid 
hydrodynamics  reduces  to  the  solution  of  the  Wave  Equation 


v2>P(?,t)  =  t) 

c  at2 


(2.5) 


in  a  linear  approximation  with  isentropic  conditions.  0 

The  constant  c  is  the  propagation  speed  and  is  defined  as 
the  partial  derivative  of  the  pressure  with  respect  to 
density  at  constant  entropy.  Both  the  pressure  and 
velocity  fields  must  satisfy  this  equation. 

By  expressing  the  field  $(?,t)  with  its  temporal 


Fourier  Transform 


P+°° 

*(r,t)  =  *(?,«) 

I  —  oo 


e-iutdco 


(2.6)  , 


the  wave  equation  reduces  to  the  Helmholtz  Equation, 


v2'V  (  r  ,  u )  +  [^]%(r,u)  =  0 


(2.7) 


For  the  remainder  of  this  work,  only  solutions  at  a  single 
frequency  u  will  be  considered;  therefore,  for  convenience, 


the  dependence  of  t  on  «  is  suppressed,  and  the  character 
-istic  wavenumber  k  replaces  the  equivalent  constant  u/c. 
This  simplifies  the  appearance  of  equation  (2.7)  to 


-  « -  *  j-  •  * 


/• 
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V2t(r)  +  k2t(r)  =  0 


(2.8)  . 


The  problems  of  acoustic  holography  either  involve  finding 
solutions  to  the  Helmholtz  Equation,  equation  (2.8),  which 
satisfy  a  particular  set  of  boundary  conditions,  or  they 
involve  finding  a  set  of  boundary  conditions  which  produce  a 
field  t  which  satisfies  equation  (2.8)  and  is  equivalent, 
over  a  surface  above  the  boundary  surface,  to  a  measured 
field. 


B .  Green ' s  Functions 


The  solution  techniques  presented  in  chapters  III-VI 
are  numerical  approximations  to  Green's  function  techniques 


for  solving  differential  equations 


The  application  of 


these  methods  to  the  solution  of  the  Helmholtz  Equation  with 
boundary  conditions  is  reviewed  in  this  section. 


B.l  Green's  Theorem  and  the  Green's  Function 


The  Green's  function  method  for  solving  the  Helmholtz 
Equation  with  boundary  conditions  involves  finding  a 
solution  to  the  inhomogeneous  form  of  equation  (2.8), 


V2G(rQ-r)  +  k2G(rQ-r)  =  -8(rQ-r) 


2.9)  , 


where  the  right-hand  side  term  is  the  delta  function,  and 


the  solution  G  is  the  Green's  function  for  a  specified 
boundary  condition.  The  utility  of  a  solution  to  equation 
(2.9)  for  solving  equation  (2.8)  is  clear  if  Green's 
theorem  is  applied  to  the  functions  G  and  t.  Green's 
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theorem  guarauntees  that 


[G(?o-?)v2*(?) 


u 


-  f (?) vaG(r0-r) ]dV  = 

0[G  {?Q-r  )  v+(  r  )  -  4*(? )  VG  (  rQ-r )  ]  4da 
S 


(2.10)  , 


where  the  right-hand  side  integration  surface  S  encloses  the 


left-hand  side  integration  volume,  and  the  surface  normal 
points  into  the  volume.  For  the  moment,  it  will  be  useful 
to  consider  f  to  be  a  solution  of  the  inhomogeneous 
Helmholtz  equation  given  by 

v2Y(3)  +  k2>f(r)  =  -f  (r)  (2.11)  . 

With  equation  (2.11)  and  (2.9),  the  left-hand  side  of 
equation  (2.10)  reduces  to 


-[k2G(r  -r)  +  v2G(r  -r) ]t(?)dV  - 


G(?o-?)f (?)dV  = 


8  (  r  - r ) f ( r ) dV  - 


G(r  -r)f(r)dV  =  +  (r„)  - 


G(?0-r)f (r)dV 


(2.12), 


so  that  the  entire  equation  simplifies  to 


*(*0)  = 


o  p 

G("?o-r)f(r)dV  +  Ofttr)  vG(rQ-r)  -  v>f  (  r )  G  (  rQ -r )  ]  °nda 

J  o 

a  (2.13). 


Equation  (2.13)  gives  an  explicit  solution  to  equation 
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i 


(2.11)  inside  the  volume  bounded  by  the  integration  surface 
S  provided  t  and  are  known  on  this  boundary,  and  the 


forcing  function  -f(?)  is  known  throughout  the  volume. 

Only  numerical  solutions  to  the  homogeneous  Helmholtz 
Equation  are  considered  in  this  work;  however,  the 
inhomogeneous  term  has  been  introduced  to  motivate  the 
existence  of  a  solution  to  equation  (2.9)  which,  along  with 
its  normal  derivative,  vanishes  over  the  integration  surface 
of  equation  (2.13)  when  this  surface  is  very  far  from  any 
volume  source  region.  It  can  be  shown  that  there  are  two 


such  solutions 


possible , 
G(?0-r) 


O  Q 

and  they  are4 

+ikR  -IkR 

:  £ -  or  “ - 

4rrR  4nR 


(2.14)  , 


where  R  =  |?o~r| .  By  requiring  the  Green's  function  to 


adhere  to  the  principle  of  causality  with  the  specified  time 
dependence  of  all  solutions  given  in  equation  (2.6),  the 
second  possibility  for  the  right-hand  side  of  equation 
(2.14)  is  ruled  out.  The  correct  form  for  the  Green's 
function  for  an  unbounded  region,  the  free-space  Green's 

function,  is  25 


3  -  e 


+  ikR 


Gf(?  -?)  = 


4  nR 


(2.15), 


and  +,  with  the  restriction  to  a  radiating  form  at  large 


distances,  has  the  solution 


*(rQ)  =  J  Gf ( f Q-r ) f ( r ) dV 


(2.16) 


It  is  shown  in  the  next  section  how  the  free-space  Green's 
function  is  of  great  utility  in  bounded  acoustic  problems 
without  volume  sources. 


.2.  The  Helmholtz  Intearals 


The  relationship  stated  in  equation  (2.13)  holds  for 
any  function  G  which  satisfies  equation  (2.9),  although  it 
is  most  useful  when  the  function  G  also  satisfies  either 
homogeneous  Dirichlet  or  Neumann  boundary  conditions.  In 
practice,  it  is  impossible  to  obtain  a  simple  closed  form 
function  satisfying  both  equation  (2.9)  and  a  homogeneous 
boundary  condition  for  anything  more  than  the  most  simple 
boundary  shape.  However,  boundary  value  problems  for 
arbitrarily  shaped  surfaces  can  be  solved  with  a  system  of 
integral  equations  formed  from  equation  (2.13)  and  the 
free-space  Green's  function. 

In  the  absence  of  volume  sources,  equation  (2.13)  with 
the  free-space  Green's  function  becomes 

+  (^0)  =  j)[t(r)  VGf  (?Q-r)  -  v*(r)Gf  (?0~r)  ]  °3da 

=  0[*o(?)9Gf(?o-?)  -  i¥N(?)Gf(?o-*)]da  (2.17) 

J  S  ^ 

where  <f  represents  the  Dirichlet  boundary  condition,  where 


^  =  ^7,  and  where  the  Neumann  boundary  condition  t N  has 

been  defined  as  -i^-j  in  anticipation  of  the  identification 

of  the  pressure  field  p  and  the  related  velocity  field 
[equation  (2.4)]  as  useful  boundary  conditions  in 
experimental  acoustics  problems.  Equation  (2.17)  is  known 

as  the  Helmholtz  Integral.29 

This  equation  alone  is  inadequate  for  solving  boundary 
value  problems,  as  it  requires  the  specification  of  both  the 
Dirichlet  and  Neumann  boundary  conditions.  This  problem  is 
readily  resolved  with  the  help  of  two  more  integral 
equations  which  relate  the  two  boundary  conditions,  and 
which  can  be  derived  from  Green's  Theorem  with  the 
Free-space  Green's  function.  However,  before  reviewing 
these  equations,  it  is  best  to  introduce  the  terminology 
used  throughout  the  remainder  of  this  work  to  distinguish 
the  region  of  the  integration  volume  in  equation  (2.12)  from 
the  rest  of  space. 

Implicit  in  equation  (2.17)  is  the  fact  that  the  point 

?  is  within  the  integration  volume  of  equation  (2.12). 

Only  radiating  problems  are  considered  here,  and  for  all  the 
boundary  geometries  considered,  a  portion  of  the  boundary 
surface  is  off  at  infinity.  The  portion  of  at  infinity  is 
an  abstraction,  of  course,  and  so  the  main  boundary  surface 
of  interest  divides  space  into  a  generally  finite  bounded 


interior  region  and  an  infinite  exterior  region.  The 
integration  volume  in  equation  (2,12)  is  taken  as  this 
exterior  region  for  radiation  problems.  To  distinguish 
between  points  located  in  the  two  regions  or  on  the 
separating  surface,  points  in  the  interior  are  subscripted 
with  I's,  points  on  the  surface  are  subcripted  with  s's,  and 
points  in  the  exterior  field  are  subscripted  with  f’s. 

With  this  terminology  established,  it  is  possible  to 
introduce  the  Interior  Helmholtz  Integral  Equation.  If  the 


point  rg  is  an  interior  point  rj,  and  hence,  not  included  in 

the  integration  volume  of  equation  (2.12),  the  left-hand 
side  of  equation  (2.10)  is  zero  in  the  absence  of  volume 
sources.  Instead  of  arriving  at  equation  (2.17)  with  the 
introduction  of  the  free-space  Green's  function  in  the 
remaining  right-hand  side  of  equation  (2.10),  the  Interior 
Helmholtz  Integral  Equation  is  obtained: 


0  =  4) 


[V*s)§£f<V*s) 


-  iH*N(rg)Gf(rI-rs)  ]da 


(2.18] 


This  equation  relates  the  Dirichlet  and  Neumann  boundary 
conditions . 


By  allowing  an  exterior  field  point  or  an  interior 


point  Irj  to  approach  the  boundary  surface  it  is  possible  to 


obtain  another  integral  relationship  between  the  two 
boundary  conditions  from  the  resulting  limit  form  of 
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equation  (2.17)  or  (2.18)  respectively.  Provided  the  normal 

direction  is  defined  at  a  surface  point  rs,  and  +  is 

reasonably  well  behaved,  it  is  readily  shown  that,  in  the 
limit,  these  equations  both  reduce  to 


V*;>  =  2J  V?s>§gf  <?s-?s>da  -  210  tN(?s,Gf  (?s"?s)da 
S-e2  S 

(2.19) , 

where  the  first  integral  excludes  an  infinitesimal  area  e2 


around  the  surface  point  .  This  Integral  equation, 

relating  the  boundary  conditions,  is  known  as  the  Surface 
Helmholtz  Integral  Equation.  The  use  of  equation  (2.18)  and 
(2.19)  together  to  uniquely  relate  the  Dirichlet  and  Neumann 
boundary  conditions  for  a  general  surface  is  discussed  in 

O  Q 

Schenck's  paper  and  in  chapter  V. 


B . 3  Rayleigh's  Integrals  and  Fourier  Transforms 


For  the  special  case  of  an  infinite  half-space  enclosed 
by  a  plane  boundary  topped  with  a  hemisphere  of  infinite 
radius,  the  solution  of  the  Helmholtz  Equation  is  much 
easier  to  calculate.  The  symmetry  of  an  infinite  plane 
makes  it  possible  to  construct  Green's  functions  which 
satisfy  equation  (2.9)  in  the  half-space  above  the  plane, 
and  which  also  satisfy  either  the  homogeneous  Neumann  or 


I 


Dirichlet  boundary  conditions  on  the  plane.  These  Green's 
functions  can  be  used  to  generate  radiating  solutions  to  the 
Helmholtz  equation  when  all  sources  are  below  the  plane, 
since  such  solutions  must  diminish  to  the  point  where 
integration  over  the  hemisphere  can  be  neglected.  On  a 
plane  boundary  in  a  Cartesian  coordinate  system  where 
H*  (x,y)  or  t  (x,y)  on  the  plane  z  =  0  is  the  boundary 

condition,  and  where  there  are  no  volume  sources,  the 
solution  to  the  Helmholz  Equation,  as  given  by  equation 
(2.13),  reduces  to  a  two-dimensional  convolution  integral: 


+<»  +~ 


*(x,y,z)  = 


, N (x‘ ,y ' )Gd , N (x~x' »Y-y ' rZ)dx'dy 1 


2 . 20) 


The  Green  functions,  G„  and  G  .  are  known: 
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G  fx.y.z)  =  Stlrlkg)6 


ikR 


2nR' 


(2.21) 


and 


_ieikR 

G  (x,y,z)  =  - 

N  2  ttR 


(2.22 


where  R2  =  x2+y2+z2.  Historically,  equation  (2.20)  for  the 
two  boundary  conditions  (D  or  N)  are  referred  to  as  the 

3  1 

first  and  second  Rayleigh  integrals. 

For  numerical  processing,  it  is  well  known  that  a 
convolution  is  more  readily  evaluated  by  using  Fourier 
32 


Transforms . 


?  will  denote  the  continuous,  infinite 


two-dimensional  Fourier  transform  of  a  function  f(x,y)  as 
defined  by: 

F(kx,ky)  =  f[f(x.y)] 

-i (k  x'+k  y ' ) 

f ( x ' , y 1 ) e  x  Y  dx'dy'  (2.23), 

and  the  inverse  operation  will  be  indicated  as  &'1.  The 
domain  of  (kx,ky)  will  be  referred  to  as  k-space ,  while  that 

of  (x,y)  will  be  referred  to  as  real-space. 

Applying  the  convolution  theorem  to  equation  (2.20) 
yields : 

4,(x,y  ,2)  =  3r"1[t0(N(kx,ky)GDN(kx,ky,z)]  (2.24). 
The  Fourier  transforms  of  the  Green's  functions  G  and  G  in 

D  N 

equations  (2.21)  and  (2.22)  may  be  found  analytically:34 

I"  exp  iz  J  k2-k2-k2  when  k2+ky<k2 

G0(kx*ky<z>  =  r  7 _ i  (2.25) 

[  exp  -z  J  k2+k2-k2  when  k2+ky>k2 

when  k2+k2gk2 

(2.26)  . 

when  k2+k2>k2 

/  k2+k2-k2 


The  circle  k2+k2=k2  is  referred  to  as  the  radiation  circle. 
For  points  (k  ,k  )  inside  the  radiation  circle,  G 


represents  the  z-direction  phase  change  of  propagating  plane 

waves,  while  for  points  (k„,k,,)  outside  the  radiation  circle 

*  y 

Gq  n  represents  the  rapid  exponential  decay  of  evanescent 
waves.  By  defining  a  complex  function  ,  kz(kx,ky)= 

J  k2-kx~ky  ,  the  transformed  Green  functions  may  be  written 
in  the  condensed  forms  G  =  exp(ik_z)  and  G  =  exp( ik„z ) /k_ . 

0  Z  N  Z  Z 

For  holography  or  inverse  scattering  problems,  the  use 
of  equation  (2.24)  is  essential,  because  in  this  way  the 
convolution  integral  (2.20)  may  be  quickly  inverted;  that 
is,  the  field  f(x,y,zii)  may  be  deconvolved  to  yield  the 

source  field  t  (x,y).  Assuming  that  t(x,y,z  )  is  known, 
we  have,  from  equation  (2.24), 

-  5-iimx.vz»)a;;»(1v|vz)1  l2-27)' 

where  G' 1 u ( k„ , k„ , zu ) ,  the  reciprocal  of  the  formula  in 

equations  (2.25)  or  (2.26),  provides  the  reconstruction  of 
both  propagating  and  evanescent  waves.  Once  the  source 
field  kt,o  N  has  been  determined,  all  other  properties  of  the 

field  may  be  calculated.^  It  should  be  noted  that  as  far 
as  the  processing  theory  is  concerned,  4*  (x,y)  need  not 

correspond  to  a  physical  source  surface.  In  fact  equation 
(2.20)  may  be  generalized  to  apply  to  the  field  and  its 
derivative  between  any  two  planes,  one  at  z  and  the  other  at 
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z  lying  between  zero  (the  actual  source)  and  z.  The  more 
general  equations  are: 


t(x, y,z)= 


4,(x'  ,y ' 


zQ ) Gq (x-x ' , y-y ' , z-zo)dx'dy' 


2.28), 


and 


( x ,  y ,  z  )  = 


+  OO  +  OO 


-i^|(x' ,y ' ,zQ )Gn(x-x' ,  y-y ' , z-zQ ) dx 1 dy ' 

(2.29). 


Equations  (2.28)  and  (2.29)  are  the  basis  of  planar 
acoustic  holography. 


C.  Definition  of  the  problems 

i 


The  first  two  sections  of  this  chapter  review  the 
fundamental  equations  of  Nearfield  Acoustic  Holography. 

This  section  defines  the  problems  of  Nearfield 
Acoustic  Holography  treated  in  the  remaining  chapters. 

C.l.  Planar  Holography 

The  evaluation  or  inversion  of  two  equations  (2.28)  and 
(2.29)  represent  four  operations  referred  to  as: 

a)  propagation  of  a  wavefront,  where  ^(x.y.z)  is 
determined  from  +  at  z 


b)  transformation  and  propagation,  where  +(x,y,z)  is 
determined  from  dv/dz  at  zQ 

c)  reconstruction,  where  f(x,y,zo)  is  determined 
from  +(x,y,z) 

d)  reconstruction  and  inverse  transformation,  where 
^{x,y,zQ)  is  determined  from  +(x,y,z). 

Since  by  definition  z>z  >0  and  zero  is  the  actual  source 

plane,  the  reconstruction  processes  may  be  referred  to 
as  the  inverse  propagation  of  a  wavefront  back  toward  the 
source.  When  z  =  zq  in  processes  b  and  d,  the  process  will 

be  referred  to  as  a  transformation. 

In  chapter  III,  the  approximations  and  assumptions 
necessary  to  reduce  the  infinite  and  continuous  convolution 
integrals  encountered  in  problems  a  and  b  to  a  finite  and 
discrete  form,  suitable  for  high  speed  numerical  processing 
are  illuminated  theoretically  and  tested  numerically.  To 
evaluate  the  convolution  integrals  two  assumptions  are  made 
first,  the  boundary  field  may  be  replaced  with  a  patch-wise 
constant  field  for  reasonably  small  patches;  and  second,  thi 
field  is  negligible  outside  of  a  finite  region.  With  these 
two  assumptions,  the  problem  reduces  to  one  of  representing 


the  Green's  functions  and  G  .  Six  methods  of  sampling  or 

D  N 

representing  the  Green's  functions  are  developed,  and  the 
effectiveness  of  these  methods  are  compared  theoretically 
and  numerically. 

In  chapter  IV,  a  numerical  inversion  method,  based  on 
equation  (2.27),  to  solve  problems  c  and  d  is  examined 
theoretically.  The  suitability  of  the  reciprocals  of  the 
six  Green's  function  sampling  forms  developed  in  chapter  III 
for  use  as  inverse  propagators  with  this  method  is  examined 
numerically.  An  iterative  deconvolution  method  is  also 
presented  in  chapter  IV. 

C . 2  Non-planar  Holography 


In  the  most  general  holographic  experiment,  neither  the 
source  boundary  nor  the  surface  over  which  field  data  is 
gathered  is  flat  in  cartesian  coordinates  or  any  other 
separable  coordinate  system.  Representing  these  surfaces 
themselves  is  an  important  problem  which  is  not  encountered 
in  planar  holography.  Chapter  V  presents  suitable  methods 
for  modelling  odd-shaped  surfaces  numerically  and  for 
relating  the  Dirichlet  and  Neumann  boundary  conditions. 
Numerical  approximation  of  the  three  Helmholtz  Integrals 
form  the  basis  of  chapter  V  and  VI. 


Chapter  VI  presents  a  computational  method  for 
determing  the  boundary  surface  conditions  of  an  odd-shaped 
radiator  from  field  measurements  above  the  surface.  The 
results  of  simulated  holograpy  experiments  which  test  this 
method  are  also  presented.  The  exact  nature  of  these 
simulations  is  fully  explained  in  chapter  VI. 
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Chapter  III. 

THE  PROPAGATION  OF  PLANAR  HOLOGRAMS 

In  Nearfield  Acoustic  Holography,  a  set  of  values, 
either  experimentally  measured  or  theoretically  generated, 
of  a  field  or  its  normal  derivative,  at  a  number  of  points 
evenly  distributed  over  a  planar,  rectangular  grid  is 
referred  to  as  a  planar  hologram.  As  representations  of 
infinite  plane  boundary  conditions,  these  holograms,  with 
the  aid  of  Rayleigh's  Integrals,  equations  (2.28)  and 
(2.29),  can  be  used  to  calculate  a  representation  of  the 
entire  three-dimensional  field  above  the  hologram  plane.  In 
this  chapter,  a  computational  method  for  determining  this 
3-D  field  is  developed,  and  test  results  of  this  method  are 
reported . 


A.  Reduction  to  finite  and  discrete  operations 

The  infinite,  continuous  integral  operations  reviewed 
in  the  last  chapter  must  be  reduced  to  finite  and  discrete 
approximations  to  allow  processing  with  current  digital 
computing  equipment.  The  validity  of  a  solution  obtained 
with  such  an  approximation  depends  directly  on  the  validity 
of  the  assumptions  and  method  used  to  reduce  the  infinite, 


continuous  relationships.  It  is  the  purpose  of  this  section 
to  clearly  delineate  the  methods  and  assumptions  used  to 
reduce  the  propagation  of  planar  holograms  to  finite  and 
discrete  operations. 


A . 1  Rayleigh's  Integrals 


The  limits  required  to  numerically  evaluate  the 
formulas  in  equation  (2.28)  or  (2.29)  with  finite,  discrete 
operations  may  be  imposed  by  restrictions  on  experimental 
data  acquisition  or  restrictions  on  computation  time  and 
capacity.  In  either  case  several  basic  assumptions  must  be 
made.  The  first  is  that  the  sources  of  the  wavefield  f  are 
such  that  the  boundary  data  H*  (x,y)  or  (x,y)  is  negligible 

outside  of  some  finite  domain  in  real  space.  Denoting  this 
domain  by  -L/2<x<L/2  and  -L/2<y<L/2,  equation  (2.20) 
[likewise  (2.28)  and  (2.29)]  becomes: 


f  ( x ,  y ,  z ) 


,+L  p + T* 
2  2 

-L  -L 
2  2 


+0 ( N (x1 ,y 1 )Gdn(x-x' , y-y 1 ,z)dx'dy'  (3.1) 


The  second  assumption  is  that  the  wavefield  t  can  be 
represented  sufficiently  well  with  a  discrete,  as  well  as 
finite,  set  of  numbers.  This  set  of  numbers  may  be  a  data 
set  from  experimental  measurements  at  lattice  points  in 
real-space  or  coefficients  for  a  superposition  of  basis 


j 


functions.  For  this  model,  the  simplest  case  is  chosen  in 
which  it  is  assumed  that  the  wavefield  t  is  not 
significantly  altered  if  the  boundary  field  ,  is  replaced 

with  a  piecewise  constant  field.  That  is,  the  L  x  L  real 

space  domain  is  divided  into  N2  patches  of  size  (L/N)  x 
(L/N),  labeled  with  integers  l,m  =  0,1,..., N-l,  and  the 
discrete  set  of  data  t  (i,m)  is  assumed  to  be  the  average 

D  /  N 

of  the  actual  boundary  field  over  the  patch.  The  field 
radiated  from  the  piecewise  constant  boundary  field  can  be 
assumed  to  represent  t  if  N  is  sufficiently  large. 

For  the  piecewise  constant  boundary  field,  equation 
(3.1)  becomes : 


N-l  N-l 

pX/+A/a  „ym+A/2 

+  ( x ,  y ,  z  ] 

i-  2  2  +D  NU,m) 

1=0  m=0  D,N 

GD,N(X-X' ,y-y' 
j  x;-a/2‘*  ym-A/2 

, z ) dx ' dy ' 

(3.2)  , 

with 

X 

ii 

( /+1/2-N/2) A 

(3.3) 

3 

II 

(m+l/2-N/2 ) A 

(3.4) 

where  A  =  L/N. 

Finally  it  is  assumed  that  it  is  sufficient  to  examine 
the  field  +  in  any  plane  z  at  discrete  points  (x  ,y  ) 

Jr  4 

defined  as  in  equations  (3.3)  and  (3.4)  for  integers  p,  q  = 
0....N  -  1.  By  defining  variables  u  =  x  -x',  v  =  y  -y', 

r  4 


equation  (3.2)  becomes: 


(3. 


„x  j  +  A/a rym+A/2 


X;~A/2 


Gd,n(xp'x' 'Yq-V .z)dx'dy' 

^m-A/a 


(  p-  1  +  1  f  2.  )  A 


( q-m+t /2  )  A 


G  M (u, v, z)dudv 

U  /  N 


d (p-i-i/2)A“ (q-m-i/2)A 
=  5D(N{p-i,q-m,z) 

with  integers  1,  m,  p,  q  =  Letting  >f(p,q,z) 

+ ( x_ , y_ , z ) ,  the  expression  for  the  radiated  field  becomes: 

r  SI 


N-l  N-l 

H-(p,q,z)  =  2  X  t.  (l,m)3  (p-J,  q-m ,  z )  (3.6), 

1*0  m=0  N  N 


which  is  the  discrete,  finite  version  of  the  Rayleigh 
integral  (2.20) . 


A. 2  Discrete  convolution 


As  with  the  infinite,  continuous  convolution,  the 
finite,  discrete  convolution  in  equation  (3.6)  is  most 
readily  evaluated  using  the  convolution  theorem 
and  Fourier  transforms  (in  particular,  the  Fast  Fourier 
Transform  algorithm,  FFT) .  Although  discussions  of  the 
discrete,  finite  convolution  theorem  may  be  found  in 
textbooks  on  signal  processing,  most  treatments  assume  that 
both  arrays  to  be  convolved  are  either  periodic  or  zero  for 
indices  outside  the  range  0,...N-1.  However,  in  the 


discrete  convolution  in  equation  (3.6)  one  of  the  arrays, 

G  (p-l,q-m,z)  [defined  in  equation  (3.5)]  may  be  evaluated 
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for  all  integers  1,  m,  p,  q,  because  the  Green's  function 
Gn.  ^(u.v.z)  is  known  analytically  over  an  infinite  domain. 

D  t  N 

Because  of  this  exceptional  feature,  the  discrete 
convolution  theorem  must  be  rederived;  the  details  of  this 
procedure  are  presented  in  Appendix  A.  Briefly,  the  results 
are  as  follows: 

The  sequence  4*  (  J  ,m)  which  is  defined  only  for 

0  i  H 

integers  (/,m)  in  (0,N-1)  must  be  extended  over  a  (2N)  x 
(2N)  domain  by  adding  zeroes  This  new  sequence  is  defined 


with : 


f  *  (i,m)  if  0</<N  and  0<m<N 

(i.m)  =  {  N 

0 '  l  0  otherwise 


(3.7) 


The  discrete  Fourier  transform  DFT(f)  of  an  array  f(/,m) 
with  (7,m)  in  (0.2N-1)  is  defined  (for  ji ,  v~0 , 1 ,  .  .  .  2N-1 )  as: 


2N-1  2N-1  -iir(/>i+mv) 

DFT(f )  =  2  Z  f ( /,m)  e  “TT 

7=0  m=0 


(3.8)  , 


and  the  inverse  is  indicated  as 

,  2N-1  2N-1  +irr (  iu+mu ) 

IDFT(F)  .  =  — =—  2  I  F(M,v)  e  (3.9). 

4N2  ji=0  u=0 

As  shown  in  Appendix  A  the  finite  discrete  convolution  (3.6) 


can  be  calculated  with: 


where 


*{p,q,z)  =  IDFT(DFT[t;N]  goN(z)) 


gn  (z)  =  DFT [ G 1  (z) ] 

3  0  ,  N  )IV  1  0  ,  N  '  '  J  )ll> 


(3.10) 


I 


(3.11) 


'S„  (1, m,z)  if  0£i<  N  and  0£m<  N 

G_  (7-2N,m.z)  if  N</<2N  and  Ogm<  N 

5*  ( i ,m, z)  =  _D'N  (3.12) 

0,N  G  ( 1.  m-2N , z )  if  0</<  N  and  N£m<2N 
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3  ( 7-2N,m-2N,z)  if  N</<2N  and  N<m<2N 

0  ,  N 

A. 3.  The  Green's  Functions 

Equations  (3.10)  and  (3.11)  form  the  computational 
basis  for  wavefront  propagation  in  which  <v(p,q,z)  is 
determined  in  any  plane  z>0  given  the  data  h*  ,(l,m),  (;, m) 

=  0....N-1.  The  data  manipulation  simply  involves  adding 
zeroes  to  the  data  sets  to  form  an  extended  (2N)  x  (2N) 
domain,  performing  a  two-dimensional  DFT  on  the  extended 
data  set,  multiplying  by  a  "propagator"  g„  (z),  and  then 

performing  an  IDFT.  The  computational  difficulty  lies  in 
evaluating  the  propagator  array  gQ  N(z)>ll;  which  must  be 

accomplished  for  each  value  of  z  and  k  (if  the  radiation 
involves  different  wavelengths  A).  If  equation  (3.10)  is 
used,  then  it  is  necessary  to  evaluate  the  integral  in 

equation  (3.5)  which  defines  .  Since  the  integration 

cannot  be  performed  analytically,  it  must  be-evaluated  with 
some  approximation.  Two  easy  ways  of  doing  this  are: 


Replace  G  <u,v,z)  in  the  integral  in  equation  (3.5) 

by  its  value  at  the  center  of  the  integration  patch; 
i.e.,  G  ( u ,  v ,  z )  a:  G  [ ( p-1 ) A , ( q-m ) A ] .  The  resulting 

1/  /  N  U  /  n 

propagator  is  as  defined  by  equations  (3.11) 

and  (3.12)  with  G  replaced  with: 

0  I  N 

S^Ur.s.z)  •  ( L/N)  2G  ( r a , sa , z )  (3.13), 

where  r,s  =  - ( N-l ) , . . . 0 , . . . N.  is  the  discrete 

transform  of  the  analytic  Green's  function  sampled  at 
discrete  points  in  real-space.  The  Green's  function 
G  (rA,sA,z)  is  singular  where  r=s=z=0,  so  when 

r  &  s  -  0  and  z  is  less  than  one  t^nth  of  a  wavelength, 
an  analytic  integrated  average  of  the  Green's  function 
is  used.  This  integration  is  over  a  circular  area 

equal  in  size  to  the  sample  lattice  square  area,  A2; 
the  details  are  in  Appendix  B.  This  form  is  refered 
to  as  the  "sampled  real-space  Green's  function." 

Replace  G„  (u,v,z)  in  the  integral  of  equation  (3.5) 

with  a  polynomial  expansion  about  the  center  of  the 
patch.  A  procedure  for  accomplishing  this  is  presented 

in  Appendix  B.  The  resulting  approximation  for  is 

given  the  symbol  5^2^(r,s,z)  with  r , s=- ( N-l ) , . . , 0 , . . , N . 

D  /  N 


gU)  (z)  is  as  defined  by  equations  (3.11)  and  (3.12) 

D  t  N 

with  G _  replaced  by  G this  form  is  referred  to  as 
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the  "integrated  real-space  Green's  function." 

In  the  derivation  of  equations  (3.10)  and  (3.11)  the 
definitions  of  the  real-space  sampling  points  given  by 
equations  (3.3)  and  (3.4)  were  used.  Because  the  DFT 
requires  an  even  number  of  sampling  points,  the  definitions 
(3.3)  and  (3.4)  result  in  an  asymmetric  distribution  of 
points  in  the  sampling  of  the  Green's  function  [as  in 
equation  (3.5)];  that  is,  there  is  an  odd  point  at  the 
origin,  N  negative  values,  and  N-l  positive  values.  By 
redefining  the  real-space  sampling,  a  symmetric  distribution 
of  points  may  be  obtained;  the  Green's  function  will  then  be 
sampled  over  patches  centered  at  (rA  +  A/2,  sa  +  A/2).  This 
shifted  form  of  equation  (3.5)  is: 

„  ( r+1 ) A  „  ( s+1 ) A 

G  (r,s,z)  =  G  ( u , v , z ) dudv  (3.14), 

u  i  n  u  t  n 

"  rA  0  sA 

with  r,s  =  -N , . . . , 0 , . . ( N-l ) .  If  a  sampled  approximation  to 
G_  M  is  used  in  equation  (3.12),  then  one  obtains  the 

"shifted  sampled  real-space  Green's  function," 

When  an  integrated  approximation  to  Gd  n  is  used  in  equation 
(3.12),  one  obtains  the  "shifted  integrated  real-space 
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Green's  function,"  g^4^(z)  .  Use  of  the  unshifted  Green's 

functions  is  equivalent  to  replacing  the  field  t(x,y,z/  with 
a  patchwise  constant  field  where  each  sampling  point  is 
centered  in  a  patch;  use  of  the  shifted  form  corresponds  to 
having  the  sampling  point  at  one  corner  of  its  patch.  If 
the  same  set  of  field  values  „  (I,m)  are  used  for  the 

0  /  N 

sample  points  in  both  cases,  the  net  result  is  a  shift  in 
the  apparent  location  of  the  represented  fields.  This  shift 
should  be  considered  when  theoretical  fields  are  used,  but 
it  is  of  little  consequence  with  experimentally  measured 
fields,  because  with  an  experimentally  measured  field,  there 
is  always  an  ambiguity  of  half  a  sample  spacing  in  the  x  and 
y  directions  when  deciding  how  to  represent  the  field  with  a 
patchwise  constant  field.  Throughout  this  study,  wavefronts 
propagated  by  the  shifted  Green's  functions  are  compared  to 
theoretical  fields  shifted  -A/2  in  the  x  and  y  directions. 

In  order  to  determine  g„  ^ ( z ) ,  from  equation  (3.11), 

D  t  N 

it  is  necessary  to  perform  the  DFT  operation.  However,  in 
the  infinite  and  continuous  formulation  in  equations 
(2.24-2.26)  the  Fourier  transformed  Green's  function 

G  (k„,k,,,z)  has  a  known  (fairly  simple)  closed  form.  If 
o  t  n  a  y 

an  approximation  for  g  (z)  can  be  found  in  terms  of 

ao,N<kx'ky-z>  in  e<3uati°ns  (2.25)  and  (2.26),  then  gDN(z)^v 
could  be  calculated  without  having  to  perform  a  DFT. 


In  order  to  derive  such  an  approximation,  equation 


(3.11)  giving  go  f z),  is  expanded  using  equations  (3.5), 
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(3.8),  (3.9),  and  (3.12),  and  the  periodicity  of  the  Fourier 

coeif  f  icients  ,  exp( -injil/N]  =  exp(  -  inji(  7-2N)  /N]  is  used. 
Equation  (3.11)  can  then  be  written  as: 


N-l  N-l  -irr(ru+sui) 

go  N(z,uu=  2  2  e  ~R 

o,n  )iv  r=_N  S=_N 


( r+i / z  )  A 


„  ( s  +  t / 2  )  A 


GD) n (x , y , z ) dxdy 


r- i /2 ) A  (s-i/z)A 


(3.15) 


The  Fourier  transform  of  is  expressed  explicilty  as: 
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GD,N^kx'ky'2)  = 


+00(,+00  -i(k  x+k  y) 

e  x  Y  GD) N(x,y,z)dxdy  (3.16) 


In  the  following  steps,  the  right-hand  side  of  equation 
(3.16)  will  be  approximated  and  compared  with  equation 

(3.15),  providing  a  relationship  between  the  G  and  g 

0  #  N  D  /  N 

The  integral  of  equation  (3.16)  is  broken  into  a  sum  of 
patch  integrals  where  exp [ -i ( kxx+k^y ) ]  is  approximated  as 

constant  over  each  patch.  The  result  is  G  a:  G 1  where: 

0  ,  N  D  ,  N 


G  d , n  ^ kx ' ky '  z  ^ 

■»  oo  -in(ru+su)  <i(r+1-/2)A(i(s  +  1/  2)A 

2  2  e  T7 

r=-<»  s=— oo 


Go , N (x' Y' z>dxdY 


(r-i/z)A  (s-i/2)A 


(3.17). 


An  estimate  of  the  error  from  approximating  exp[ -i ( k„x+k„y ) ] 

x  y 

as  constant  over  each  patch  is  shown  in  Appendix  C  to  be: 
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kxA  kyA 

Go,N(kx'ky'z> 

4  sin(kxA/ 2 )  sin( kyA/ 2 ) 

(3.18) 


This  error  is  limited  if  the  2N  values  of  kx  and  ky  are 


chosen  so  as  to  minimize  the  distance  from  the  origin  in 
k-space .  This  restriction  along  with  consideration  of  the 
Fourier  coefficients  in  equation  (3.15)  motivates  the 
definitions : 


_  f  M 


if  0  <  ji  <  N 


xM 


.  I  Ni 

l  n(  u-2N)  if  N  S  ji  < 


(3.19) 


2N 


and  similarly  ky^  in  terms  of  v .  Evaluating  (kx>ky)  at 


(kx^,kyv)  and  using  equation  (3.5)  defining  Gd  n,  equation 


(3.17)  can  be  written: 


G„  (kv  ,k„  ,z) 

0  ,  N  '  Xji  y  v  ' 


N-l  N-l  -in  ( rji+su)  00  _ 

2  2  e  ~N  2  2  ( r+2mN , s+2nN , z ) 

r=-N  s=-N  m=-«>  n=-»  0  '  N 

(3.20) 


The  m=n=0  term  is  simply  g„  (z)  ,  so  that  equation  (3.20; 
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can  be  written  as: 


G  ( k„  ,  k„  ,  Z  )  s: 
o  ,  n  x  ji  y  v 


g 


0  ,  N 


N-l  N-l  -  in  (  r  ji+s  v )  °°  °° 

(z)  +  2  2  e  ~N  2  2  5n  „ ( r+2mN , s+2nN , z ) 

M  r=-N  s  =  -N  m=-<=°  n=-«> 

( m=n ) *0 


D  ,  N 


(3.21) 


One  should  note  that  if  Go  N(kx^,kyv,z)  replaced 


g  (z)  in  equation  (3.10),  only  the  m=n=0  term  of  equation 


(3.20)  [the  g 


(z)  term  in  equation  (3.21)]  will  yield 

the  result  desired  as  in  the  discrete  convolution  of 
equation  (3.6).  The  summmation  remaining  in  equation  (3.21) 
effectively  extends  the  range  of  the  summation  in  equation 
(3.6)  to  include  the  infinite  periodic  extension  of  the 
finite  input  data  field.  This  extension  may  be  viewed  as 
representing  an  infinite  set  of  fictitious  image  sources;  an 
apparent  copy  of  each  field  point  is  located  at  all  of  the 
points  ( x+2mL , y+2nL )  for  all  pairs  of  integers  (m,n).  These 
fictitious  image  sources  can  produce  large  errors  in 
numerical  calculations  of  propagation. 

If  G„  fx,y,z)  falls  off  rapidly  with  x  and  y,  then 

G„  (x,y,z)  also  falls  off  rapidly,  and  the  effects  of  image 

sources  will  be  small.  In  this  case  the  summation  in 
equation  (3.21)  may  be  ignored,  and  equation  (3.21)  then 

shows  that  G0  ,  N  ( kX;i  ■  ky„  -  z )  approximates  gD/N(z)  directly. 
This  approximation  is  called  the  "sampled  k-space  Green's 
function",  is  labeled  gj^fz)  ,  and  is  defined  by: 
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9d,h^z^>iv  GD,N^X)i,kyv,z* 


(3.22) 


On  the  other  hand  if  G„  (x,y,z)  does  not  fall  off 

D  /  N 


rapidly,  then  the  sum  in  equation  (3.21)  may  be  significant 
and  Gd  n  may  be  an  unacceptable  approximation  of  go  N> 


Inspection  of  equations  (3.16-3.21)  shows  that  this  sum 
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could  be  reduced  if  G„  (x,y,z)  in  equation  (3.16)  were 

D  i  N 

replaced  by  G  (x,y,z)  multiplied  by  a  suitable  window 

function  W(x,y) .  This  window  function  should  have  a  value 

near  unity  for  x  and  y  in  the  range  [-L,L]  and  should  fall 

off  rapidly  outside  this  range.  Such  a  window  would  reduce 

any  fictious  image  sources  without  introducing  any  new 

3  1  ~  • 

error.  A  simple  closed  form  function  G'  is  sought  such 

that : 


+  00  +00  _ 


GD,N(kx'ky'2>  = 


i(k  x+k  y) 
e  Y  W(x,y)G(x,y,z)  dxdy 

U  I  N 


(3.23), 

and  the  same  process  outlined  in  equations  (3.16-3.21)  yield 
an  expression  analogous  to  equation  (3.21): 

G "  ( k  ,  k„  ,  z )  a:  g  ( z )  + 

d,n'  Xji  y  v '  '  3o,n'  '  )IU 

N— 1  N— 1  -in(  ru+sv )  °°  °° 

2  2  e  “N  2  2  W( r+2mN, s+2nN) 5  ( r+2mN , s+2nN , z ) 

r=-N  s=-N  m=-°°  n=-«>  0 ' N 

(m=n) ^0  (3.24). 

Because  of  the  window  function,  this  summation  is  not  as 

significant  as  that  of  equation  (3.21)  regardless  of  how 

fast  G  (x,y,z)  decreases  with  increasing  x  and  y. 

Ideally  the  value  of  the  window  function  would  be  unity 
within  the  limited  range  (corresponding  to  m=n=0)  and  zero 
outside,  but  neither  equation  (3.23)  nor  its  equivalent 
convolution  integral  in  k-space  yield  a  simple  closed-form 
expression  for  such  a  window.  The  goal  here  is  to  reduce 


the  computation  time  of  g  (z) 

r  3  D  ,  N  '  '  )IV 


in  equation  (2.11) ; 


therefore  a  window  for  which  equation  (3.23)  can  be 
evaluated  analytically  is  necessary. 

The  "radial  sine  function,"  2Ja (kQr ) /kQr  [Ji(z)  is  the 

first  order  Bessel  function,  kQ=  n/L,  and  r2=  x2+  y2 )  ]  is  a 
suitable  window  function.  It  does  provide  a  simple  and 
closed  form  expression  for  G"  .  This  window  function 
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reduces  the  sum  in  equation  (3.23)  so  that  G"  (k„  ,k„  ,z) 
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more  accurately  approximates  go  N(z)  ^  than 
Go  N ( kx^ , ky v , z ) .  Using  the  convolution  theorem,  it  can  be 
shown  that  G"  for  this  window  is  approximately  an 
integrated  average  of  G„  ^  over  an  annular  region  around 
points  (k  ,k  ).  This  function  G"  is  the  integrated 
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k-space  Green's  function,  labeled  g^6^(z)  ;  which  is 

r  3  d  ,  n  '  '  y.v 

defined  as35 


nkr+ko  - 


g^  6  ^ (z ) 

3 o  ,  n  '  '  y.v 


2k„k„  o 


r~o 


k2 


k  -k 
r  o 


G  ( k ' , z ) k '  dk '  k  >k 

d  ,  n  '  r  ’  r  r  r-  c 


(3.25; 


nk 


°  aD,*(kr'z>kr  dkr 


kr=0 


where  k2  =  k2^  +  k2y.  Closed  form  evaluations  are  presented 


in  Appendix  C. 


significantly  as  z  varies  from  a  fraction  of  a  wavelength  to 
several  wavelengths.  This  particular  range  of  propagation 
distance  is  of  great  interest  in  NAH  applications.  In  this 
subsection,  plots  of  the  Dirichlet  and  Neumann  Green's 
functions  and  their  Fourier  Transforms  are  presented  for 
propagation  distances  of  a  fraction  of  a  wavelength  and 
several  wavelengths. 

B.l.  Dirichlet  boundary  conditions 

The  analytic  Fourier  transform  of  the  Dirichlet  Green's 
function,  equation  (2.25),  is  plotted  in  Figure  3.1a  for  the 
small  distance  z/X  =  0.13  and  extending  over  a  typical  range 
of  kx  with  ky=  0.  The  slowly-varying  nature  of  the  function 

allows  for  easy  sampling  with  a  uniform  spacing  of  points 
(kx,ky) . 

The  real-space  function,  equation  (2.21),  for  the  same 
distance  z/X  is  plotted  in  Figure  3.1b.  Since  the  source 
fields  considered  in  this  paper  are  assumed  to  be  negligible 
outside  a  region  of  length  L,  the  Green's  functions  in 
real-space  need  only  be  considered  in  the  range  ( -L , L ) . 
However,  in  contrast  to  the  k-space  form,  the  spiked  shape 


of  the  real-space  function  makes  it  difficult  to  sample  with 
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a  finite  number  of  points  with  regular  spacing;  it  becomes 
impossible  as  z  goes  to  zero  and  the  function  becomes 
singular . 

Figures  3.1a  and  3.1b  suggest  that  for  small  distances 
( z/X  <  0.5)  a  direct  sample  of  the  k-space  Green's  function 
will  yield  good  results.  This  expectation  is  also 

consistent  with  the  relationship  of  g^5^(z)  to  g  (z) 

developed  in  section  A. 3.  and  summarized  by  equations  (3.21) 
and  (3.22).  Since  most  of  the  area  under  the  curve  of 
Figure  3.1b  is  contained  in  a  relatively  small  area  around 
x=y=0,  the  sum  in  equation  (3.21),  which  quantifies  the 
fictitious  image  source  error,  will  be  very  small,  and  hence 

g„  ( z )  will  be  represented  well  by  g^5^.  For  small 

distances  z/X,  conditions  are  not  favorable  for  the  direct 
sampling  of  the  real-space  function  due  to  the  singularity 
atx=y=z=0.  Instead  of  a  direct  sampling  of  the 
real-space  function,  a  more  sophisticated  approximation  of 

the  integration  in  equation  (3.5)  is  needed,  such  as  g^2^(z) 
or  g£4) (z)  . 

Figures  3.1c  and  3.  Id  are  plots  of  the  analytic 
transform  and  real-space  form  of  the  Dirichlet  Green's 
function  for  z  greater  than  three  wavelengths  (z/X  =  3.1). 
The  analytic  transform  represents  the  real-space  Green's 
function  of  infinite  extent.  This  function  does  not  fall  off 


rapdily  in  real  space,  and  this  shows  up  in  k-space  as  the 
rapid  oscillations  of  Figure  3.1c  which  are  very  difficult 

to  sample  directly  as  in  g^5^(z).  The  Green's  function  for 
this  distance  might  be  represented  by  the  integrated  average 
k-space  form  g^6^(z);  however,  the  slow  variation  evident  in 

Figure  3 . Id  suggests  that  a  direct  sampling  of  the 
real-space  function  will  be  an  adequate  representation  of 
the  propagator.  Using  one  of  the  real-space  Green's 

function  forms,  g^^(z)  through  g^4^(z),  along  with 

augmenting  the  data  field  with  zeros  as  in  equation  (3.6), 
will  eliminate  the  fictitious  image  source  error. 

B.2.  Neumann  boundary  conditions 

The  analytic  transform  of  the  Neumann  Green's  function, 
equation  (2.26),  is  plotted  for  small  z/X  in  Figure  3.2a 
while  the  real  space  function,  equation  (2.22),  is  plotted 
in  Figure  3.2b.  Unlike  the  well-behaved  function  plotted  in 
Figure  3.1a  the  singularity  of  the  k-space  Neumann  Green's 

function  makes  it  unsuitable  for  direct  sampling  as  g^,5^(z). 
The  integrated  average,  gf,6^(z),  is  necessary  to  sample  the 

N 


Green's  function  from  k-space.  While  the  k-space  function 


k“  for  all 


is  singular  on  the  radiation  circle  k~  +  kn  = 

x  y 

distances  z/X,  the  real  space  function  is  singular  only  at 
the  point  x=y=z=0 .  This  singularity  presents  no  computa¬ 
tional  problem  if  either  g^^  (z)  or  g^4^(z)  is  used  to 
represent  g  (z). 

The  k-space  form  of  the  Green's  function  for  Neumann 
boundary  conditions  at  the  larger  distance  z  =  3A,  depicted 
in  Figure  3.2c,  shows  the  same  rapid  oscillations  as  the 
Dirichlet  version,  and  it  is  also  singular  on  the  radiation 

circle.  The  integrated  form  g^6^(z)  must  be  used  if  a 

k-space  Green's  function  is  used,  but  again,  the  slow 
variation  of  the  real  space  form,  Figure  3. 2d,  suggests 

using  one  of  the  real  space  forms  g^^(z)  through  g^4^(z). 

N  N 

The  eight  plots  of  the  Green  functions  presented  are 
representative  of  all  the  cases  studied.  They  are  intended 
to  serve  as  a  basis  for  understanding  the  results  of  the 
computational  experiments  presented  below. 

C.  Quantitative  comparison  of  the  sampled  Green's  functions 

The  purpose  of  this  section  is  to  present  a 
quantitative  comparison  of  the  effectiveness  of  the  six 

sampled  Green's  function  forms,  g^^(z)  through  g^6^(z)  when 


they  are  used  in  each  of  the  four  processes:  propagation, 
transformation,  reconstruction,  and  inverse  transformation. 
The  comparisons  span  a  suitable  range  of  separation  distance 
between  the  input  and  output  fields  for  each  case.  To 

compare  the  six  forms,  readily  calculated  fields  ^(x.y.z) 
and  ^(x,y,z)  of  a  tractable,  theoretical  boundary  condition 
4*  (x,y)  are  used.  These  fields  are  used  as  simulated  input 

N 

fields,  or  holograms,  and  as  absolute  reference  fields  in 
comparisons  of  the  six  output  fields.  An  output  field 

(x,y,2)  or  ^|x!y,z)  is  associated  with  each  sampled 

Green's  function  g^Jtz)  ,  i  =  1,...,6,  used  in  the  numerical 

processing  of  the  input  data. 

Since  the  processing  problems  which  have  been  discussed 
have  direct  relevance  to  Nearfield  Acoustic  Holography  it  is 
natural  to  test  and  compare  the  different  Green's  function 

fronts  g(*l(z)  in  simulated  experiments  with  an  acoustic 
pressure  field,  4»(x,y,z)  =  p(x,y,z).  The  Neumann  boundary 
condition,  and  -i^(x,y,z)  for  any  z,  are  specified  in  terms 
of  the  z-component  of  the  particle  velocity  defined  by 

vz<x’Y'z>  =  -j&r  fzP(x.Y-Z)  (3.26), 

where  pc  is  the  characteristic  acoustic  impedance. 
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C.l.  A  theoretically  tractable  boundary  condition:  the 
baffled  piston 


The  computationally  tractable  and  relatively 
nonsingular  acoustic  pressure  field  of  an  oscillating  piston 
in  an  infinite  baffle  is  used  for  the  NAH  simulations.  In 
particular,  the  chosen  source  is  a  baffled  piston  of  radius 
a  =  3.44/k  and  surface  velocity  amplitude  u  .  This  source 

radiates  into  a  medium  with  a  characteristic  impedance  pc. 

Calculating  the  theoretical  field  for  each  field  point 
requires  the  evaluation  of  a  one-dimensional  integral  with 
fixed  limits.  From  equations  (2.20-2.22  and  3.26)  the 
fields  produced  by  the  piston  are  given  by 


p° (x,y , z) 


i£ckuo 

2  TT 


n2n  pka 

e? 

J  0  “  0 


ikR 

R 


rdr  de 


(3.27) 


and 


v£(x,y,z)  =  — o 


„2n 


2tt 


„ka 


z ( 1-ikR ) 


,ikR 


rdr  do  (3.28), 


where  R2  =  (x2  +  y2 )  +  r2  -  2(x2  +  y2)1,/2  cos(0)+  z2 . 

Equations  (3.27)  and  (3.28)  can  be  reduced  to  the 


following  one  dimensional  integrals: 

.TT 


36 


p°  (x,y,z)=  f  eikR[ - 

2  n  1  Jn  [l  + 


b  cos ( 9  )  -  1 
b2  -  2b  cos(0) 


de  +  a(b)eik2j 


(3.29) 


2tt 


v°(x,y,z)=  -^o{  I  S  eikRr - b...cosl8)_-.JL 


[l  +  b2  -  2b  cos(0)J 


de  +  a(b)eikz} 
(3.30), 


where  b2  =  (x2  +  y2 ) /a2  and 


a  ( b )  = 


0  ,  b  >  1 
n  ,  b  =  1 
2rr  ,  b  <  1 


C.2.  Sampling  the  theoretical  field 

By  numerical  integrations  of  equations  (3.29)  or  (3.30) 
with  a  fixed  value  of  z,  the  fields  p°(x,y,z)  or  v°(x,y,z) 
are  evaluated  at  each  of  the  coordinate  locations  (x;<Ym)» 

defined  in  equations  (3.3)  and  (3.4),  yielding  a  64  x  64 
array.  Each  element  of  the  array  is  a  complex  number 
specifying  the  field  amplitude  and  phase  at  the  point 

(x^,ym/Z).  Elements  of  these  arrays  are  written  as  p°m  or 
v/m  with  value  of  z  taken  from  context  if  not  explicitly 


indicated . 


C.3.  Quantifying  the  error  in  the  output  fields 


In  order  to  evaluate  the  relative  accuracy  of  the 
different  Green's  functions  the  output  fields  4^^  processed 

with  9^1  are  compared  with  the  reference  fields  p^m  or 

v°m,  and  an  average  difference  and  standard  deviation  are 

determined.  All  quantities  are  normalized  by  the  largest 
magnitude  of  the  reference  field  in  the  array. 

The  difference  in  the  real  (Re)  and  imaginary  (Im) 
parts  at  each  point  are  defined  as 

]}  >  max [ |  +im | ]  <3'31> 


A)in  *  {Im[tJm)]  -  ]}  '  max  [  |  j  ]  <3'32> 


The  average  difference  is  defined  as 


—  fi^  f  r  R  ^  ^ )  t  ( i )  I 

a(  ’  "  {  mlQ  [A/m  +  *?m  ]  }  /  8192  <3-33> 


and  the  standard  derivation  is 


11  ■  {  jo  Jo  {Ki11- 


"I  1  /  2 

/8192 j 

(3.34) 


The  results  of  each  case  are  summarized  in  a  chart  such  as 
that  of  Figure  3.3a.  A  bar  is  drawn  for  each  function 


g'  ' (z) ,  i  =  1  to  6,  with  the  center  drawn  a  distance  A  from 
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the  zero  axis  while  the  top  and  bottom  are  drawn  a  distance 
+o  and  -a  from  the  center. 

To  further  illustrate  the  effectiveness  of  each  method, 
several  plots  are  presented  which  display  a  central  cross 
section  of  either  the  real  or  imaginary  part  of  the  output 

field  represented  by  crosses,  against  the  corresponding 

reference  field  p^m  or  v°m  represented  by  a  solid  line. 


C ■ 4 .  Test  results  for  wavefront  propagation 

In  this  section,  charts  and  plots  are  presented  which 
quantify  and  illustrate  the  effectiveness  of  the  sampled 
Green's  function  forms  when  used  to  propagate  a  pressure 
wavefront  and  to  transform  and  propagate  a  z-component  of 
the  particle  velocity  wavefront.  As  defined  in  chapter  II, 
these  processes  are  the  numerical  evaluations  of  equations 
(2.28)  and  (2.29)  for  fixed  z-zq  >  0.  Data  for  four 

different  values  of  (z-zq)/A  are  presented  (to  simplify 
notation  in  the  figures  (z-z  )/A  is  written  z/A). 

The  chart  in  figure  (3.3a)  graphically  displays  the 
average  differences  and  standard  deviations  of  the 

output  fields  p  ^ ( z=0 . 006 X )  from  the  field  p^„ ( z=0 . 006A ) 


Fig.  3.3  The  normalized  error  (defined  in  eguations  47  and 
48)  produced  by  each  of  the  six  Green's  function  sampling 
methods  when  used  to  propagate  the  pressure  field 
progressively  larger  z-distances,  (a),  (b),  (c),  and  (d) 

from  an  input  theoretical  baffled  piston  pressure  field. 


when  the  input  field  is  the  pressure  at  the  theoretical 
piston  source  boundary,  P5m^zo=0^‘  The  ^-space  Green's 
function,  equation  (2.25),  is  well  behaved  for  small  (z-zQ) 
and  so  requires  only  a  simple  representation;  figure  (3.3a) 
shows  that  g^5^(z-zQ)  and  g£6^(z-zQ)  do  yield  good  results 

for  very  small  propagation  distance.  The  real-space  Green's 
function,  equation  (2.21),  is  nearly  singular  for  small 

(z-zQ)  and  so  requires  a  sophisticated  represenation ;  figure 

(3.3a)  shows  that,  of  the  real-space  forms,  g^4^(z-zQ) 

yields  the  best  results  for  very  small  propagation 
distances . 

The  input  field  for  the  chart  of  figure  3.3b  is  again 

p°m(zQ=0)  while  the  six  output  fields  are  p^^ ( z=0 . 13X ) . 

At  this  distance,  all  six  functions  yield  good  results. 

For  the  chart  in  fiqure  3.3c,  the  input  field  is 

p^m(  zQ=0 . 13X )  and  the  output  fields  are  p^^(z=l.lX).  At 

this  distance  the  error  from  any  of  the  real-space  forms  is 
about  one  third  that  of  either  k-space  form.  There  is  little 
difference  between  the  errors  from  the  four  real-space 
forms.  In  general  at  yistances  (z-zQ)  >  X/2,  the  real  space 

forms,  g£*^(z)  through  g£4^{z),  provide  better  results  than 


the  k-space  forms . 


For  a  propagation  distance  of  several  wavelengths,  the 
chart  in  figure  3.3d  indicates  that  all  of  the  real-space 

Green's  function  sampling  forms  g^)(z)  through  g^4^(z) 

yield  similar,  good  results.  The  relative  error  of  the  the 
k-space  forms  grows  quickly  with  an  increase  in  propagation 

distance  with  g^®^  giving  the  largest  error.  This  large 

error  is  due  to  unattenuated  fictitious  image  sources. 

The  charts  of  figure  3.3  suggest  that  the  shifted, 

integrated  real-space  Green's  function  g^4^(z)  is  the  best 

choice  for  a  single  function  which  will  yield  good  results 
over  the  range  of  (z-zQ)  considered  here.  The  plots  of 

figures  3.4  and  3.5  further  illustrate  the  deficiencies  of 
the  other  functions  in  comparison  with  the  merits  inherent 

in  g£4^(z).  In  figure  3.4a  a  central  cross  section  of  the 

imaginary  part  of  the  output  field  p^^ ( z=0 . 006A )  is  plotted 

with  the  imaginary  part  of  p°m ( z=0 . 006X )  while  the  same  is 

done  for  p ^4 ^ { z=0 . 006X )  in  figure  3.4b.  In  both  cases  the 
propagation  distance  (z-zQ)  was  0.006X.  Figure  3.4a  shows 

that  underestimates  the  imaginary  part  of  the  field 

while  figure  3.4b  shows  that  p^4^  is  a  much  better 
approximation.  As  a  first  order  approximation  of  G  (z), 
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Fig.  3.5  The  real  part  of  the  pressure  3.1*  from  the  plane 
of  a  theoretical  baffled  piston  source  generated  from  an 
input  pressure  field  sample  taken  0.13  $L  from  the  source 

using  (a),  plus  signs.  gl5)(z)  and  (b)  plus  signs  g£4)(z). 

Theory  is  the  solid  line. 


equation  (23),  g^  1 (z)  apparently  cannot  represent  the 

behavior  of  the  nearly-singular  real-space  Green's  function 
at  very  small  distances.  The  higher  order  approximation 

g£4^(z)  does  appear  to  represent  the  effect  of  this  behavior 

well . 

For  propagating  over  distances  of  a  wavelength  or  more, 
the  effects  of  image  sources  should  degrade  the  performance 

of  the  k-space  forms  g^5^(z)  and  g^6^(z).  Figure  3.5a 
clearly  shows  an  interference  effect  from  these  image 
sources  in  the  plot  of  p^^(z=3.2A)  against  p°m(z=3.2A) 
where  (z-zq)  was  3.1*.  The  integrated  real-space  form  does 
not  suffer  from  any  fictitious  image  source  error  as  can  be 
seen  in  figure  3.5b  where  p^(z=3.2*)  is  compared  with 

P;m<  z=3.2*),  the  same  conditions  as  in  figure  3.5a. 

The  charts  and  plots  of  figures  3.3,  3.4  and  3.5 
compare  the  utility  of  the  six  Green's  function  sampling 
methods  when  used  in  the  numerical  evaluation  of  equation 
(2.28).  The  charts  and  plots  of  figures  3.6,  3.7,  and  3.8 
do  the  same  type  of  comparison  for  numerical  evaluations  of 
equation  (2.29).  The  chart  in  figure  3.6a  displays  the 

average  differences  and  standard  deviations  of  the 

output  fields  p^.^(z=0)  compared  to  the  reference  field 
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Fig.  3.6  The  normalized  error  produced  by  each  of  the  six 
Green's  function  sampling  methods  when  used  to  propagate  the 
pressure  field  progressively  larger  distances,  (a),  (b) , 

(c),  and  (d),  from  an  input  theoretical  z-component  of  the 
particle  velocity  boundary  field. 


P2m^z=0^‘  T^e  ^nPut  field  is  the  normal  component  of  the 
surface  velocity  at  the  theoretical  piston  source  boundary, 
v5m(zo=0) 1  Functions  g(4^(z)  and  g^6^(z)  yield  the  best 

results  for  this  transformation  from  v/m(z0=°)  to  p^^(z=0). 
Unlike  the  functions  g£^(z),  there  appears  to  be  no 

value  of  (z-zq)  where  all  g^^(z)  produce  similar  results; 
this  can  be  seen  in  figures  3.6a  through  3.6d.  The  function 
g^4^(z)  yields  as  good  results  as  any  function  over  the 
range  of  distances  (z-zQ)  considered.  The  real-space  forms 

gi2^  (z)  and  g^3^(z)  are  only  slightly  inferior  to  g|,4^(z) 

N  N  N 

for  use  when  the  distance  (z-zQ)  is  small.  The  real  space 

form  gi^(z)  starts  out  the  worst  of  all  six  forms  and 

steadily  improves  with  propagation  distance  while  the 
k-space  forms  decrease  in  performance  with  increasing 
propagation  distance. 

Figures  3.7a  and  3.7b  compare  with  p^m(z=0)  the  fields 
p^^(z=0)  and  p^4^{z=0)  respectively  as  produced  with 
g^5^(z=0)  and  g^4^(z=0).  Against  p°m(z=3.2A),  the  fields 
p<®>(z-3.2A)  and  p^4  ^  (  z=3 . 2A )  ,  as  produced  with  g^5^(z=3.2*) 
and  g^4 ^ ( z=3 . 2X ) ,  are  compared  in  figures  8a  and  8b 
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Fig.  3.7  The  real  part  of  the  pressuce  in  the  plane  of  a 
baffled  piston  source  along  a  center  row  of  the  output 
generated  field  array  produced  using  (a),  plus  signs, 

g^5^ (z)  and  (b) ,  plus  signs,  g^4^(z)  with  the  z-component  of 

N  N 

the  particle  velocity  in  the  same  plane  as  input  to  the 

?rocessing  routines.  The  solid  lines  represent  the 
heoretical  value. 


displacement  (x/X) 


Fig.  3.8  The  real  part  of  the  pressure  field  3.2  A  from 
the  input  baffled  piston  as  produced  starting  with  the 
z-component  of  velocity  at  the  source  plane  using  (a),  plus 

signs,  g^^(z)  and  (b),  plus  signs,  g^4^(z).  The  solid 

lines  plot  theory. 
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respectively.  In  both  figures  3.7a  and  3.8a,  the  effects  of 
fictitious  images  introduced  by  ;(z)  are  seen.  By 
contrast  both  3.7b  and  3.8b  show  the  good  agreement  obtained 
with  g( 4) (z) . 

N 


D.  Conclusion  and  summary 

Over  the  full  range  of  distances  (z-zq)  considered  and 

for  both  propagation  problems  and  problems  of  transformation 
and  propagation,  the  shifted  integrated  real-space  Green's 

functions  g[4l,{z)  appear  to  yield  consistently  good  results. 

This  consistency  makes  this  form  ideal  for  computation. 

Section  III  developed  the  theory  for  approximating  the 
infinite  and  continuous  integrals  with  finite  and  discrete 
forms.  The  approximation  was  based  on  two  main  assumptions: 
First,  the  boundary  field  can  be  neglected  outside  of  a 
finite  region,  and  second,  this  field  can  be  adequately 
represented  by  a  patchwise  constant  field  for  reasonably 
small  patches.  If  these  two  assumptions  are  satisfied,  it 
was  shown  that  the  problem  was  to  represent  the  known 
Green's  function  in  a  discrete  form,  and  six  methods  for 
accomplishing  this  were  developed.  A  qualitative 
examination  of  the  actual  form  of  the  Green's  functions  was 


given  in  section  IV.  This  was  done  to  identify  the 


difficulties  in  representing  the  Green's  functions  and  to 
identify  how  these  difficulties  might  be  overcome  by  one  or 
more  of  the  six  methods.  The  results  of  testing  presented 
in  section  V.E  suggest  the  use  of  the  shifted  integrated 

real  space  Green's  function,  g^4](z),  for  acoustic 

D  i  N 

propagation  problems.  The  results  show  this  form  produced 
minimum  of  error  for  propagation  distances  from  zero  to 
several  wavelengths . 
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Chapter  IV 


SOURCE  FIELD  RECONSTRUCTION  FROM  PLANAR  HOLOGRAMS 


In  this  chapter,  the  inverse  of  the  problem  considered 
in  the  last  chapter  Is  examined.  Here  it  is  assumed  that 
a  known  planar  hologram  H'(p,q,z)  accurately  represents  the 

field  in  the  hologram  plane,  and  that  the  field  between  the 
hologram  plane  and  the  source  plane,  z=0 ,  is  sought. 

In  the  first  section  of  this  chapter,  the  numerical 
reconstruction  of  source  fields  is  developed  as  a  direct 
approximation  to  the  continuous  formulation  of  equation 
(2.27).  This  is  essentially  the  technique  suggested  in  the 

original  Nearfield  Acoustic  Holography  paper1  with  minor 
improvement . 

The  second  section  presents  test  results  illustrating 
the  suitability  of  the  reciprocals  of  the  six  Green's 
function  sampling  methods,  developed  in  chapter  III,  as 
inverse  propagators  with  the  method  of  section  A. 

An  important  original  contribution  of  the  first  section 
of  this  chapter  the  definition  of  error  sources  in  the 
numerical  method  there  described.  The  identification  of 
these  errors  motivate  the  search  for  alternative 
deconvolution  methods,  and  in  the  last  section,  an 
alternative  numerical  source  field  reconstruction  method  is 
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presented.  This  method  is  based  on  the  theory  developed  in 
chapter  III  .  That  theory  led  to  equation  (3.6)  which  can  be 
written  in  matrix  form: 


N2-l 

+(r, z  )  =  £ 

H  s=0 


K,,.<z„)]r8VK<s>  f4-1 


where  the  N2  unique  index  pairs  (p,q)  and  (/,m)  of  equation 
(3.6)  are  themselves  indexed  respectively  by  the  indicies  r 

and  s,  and  where  the  elements  of  the  N2  x  N2  matrix 

[G(d,n)  (ZH>jrS  are  SD,N(P'I'q-|n'ZH1-  The  SOUrCe  field 
reconstruction  problem  would  be  solved  by  inverting  this 
matrix,  but  inverting  this  matrix  for  large  N  would  require 
excessive  computational  capabilities  and  will  preclude  easy 
high  speed  data  processing.  The  stability  of  such  an 
inversion  is  also  suspect.  However,  iterative  methods  for 
solving  this  equation  without  inverting  the  matrix  are  quite 
feasable.  Such  a  method  for  solving  equation  (4.1)  is 
presented  in  section  C. 


A.  Deconvolution  with  the  Discrete  Fourier  Transform  and 
Smoothing 


For  the  case  of  infinite  and  continuous  fields,  the 


convolution  theorem  provides  an  exact  formula  for  the 
deconvolution  or  reconstruction  of  +  (x,y)  from  t(x,y,z  ) 
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as  in  equation  (2.27).  It  would  be  tempting  to  try  the  same 
technique  with  the  discrete  convolution  given  in  equation 
(3.10).  One  could  write 

*;iN(/.a)  =  IDFT[DFT{f(p,q,zH)}/g0(N(zH)]  (4.2). 

However,  equation  (3.10)  can  be  derived  (without  approxima¬ 
tions)  only  if  the  DFT  is  defined  for  a  2N  x  2N  domain. 
Equation  (3.10)  was  found  by  assuming  that  +'  „  (;,m)  was 

zero  at  the  extra  points  in  the  extended  domain.  In  order 
to  evaluate  equation  (4.2),  the  (in  general  non-zero)  values 
of  t(p,q,z  )  must  be  known  at  the  extra  points  in  the 

extended  domain.  However,  in  holography  the  hologram  data 
is  defined  to  be  the  N  x  N  array  t(p,q,z  ).  Although  N  can 

be  redefined  so  that  the  hologram  data  covers  the  (2N)  x 
( 2N )  domain,  it  is  almost  certain  that  the  data  measured 
outside  the  N  x  N  domain  do  not  correspond  to  the  ideal 
values  compatible  with  the  derivation  of  equation  (4.2). 
Furthermore,  one  must  consider  the  possibility  that  for  some 
(ji,u)  the  array  element  g„  „  ( z )  may  be  zero;  in  the  ideal 

case  DFT{* (p,q,zu) } , ,  will  also  be  zero  and  the  quotient  in 

equation  (4.1)  will  be  indeterminate.  If  t(p,q,zu)  is  the 

data  measured  over  the  2N  x  2N  domain  it  is  very  unlikely 
that  DFT(t(p,q,z  ) }  will  be  exactly  zero,  so  that  its  use 

in  equation  (4.1)  may  lead  to  serious  computational  errors. 
Thus  for  measured  hologram  data,  equation  (4.1)  cannot  be 


Another  approach  to  the  inverse  scattering  problem  is 
to  assume  that  both  +  (x,y)  and  t(x,y,zu)  are  negligible 

for  sufficiently  large  x  and  y.  This  is  especially  valid  in 
Nearfield  Acoustic  Holography,  where  the  hologram  data  is 
measured  as  closely  as  possible  to  the  sources,  which  are 
assumed  finite  in  size.  The  measured  hologram  data  is 
assumed  to  comprise  the  N  x  N  array  t(p,q,zu) ,  for  p,q  = 

0,1,... N-l.  As  already  mentioned,  equation  (10)  provides  an 
exact  deconvolution  formula  for  infinite,  continuous  fields. 
This  equation  can  be  employed  by  approximating  the  infinite 
continuous  Fourier  Transform  of  the  hologram  field 
+(x,y,z  ).  With  the  assumed  finite  extent  of  4»(x,y,z  )  the 

DFT  of  the  hologram  data  f(p,q,zH)  is  a  logical  choice  for 

such  an  approximation,  but  the  DFT  of  a  sampled  field  will 
differ  from  the  continuous  field  transform.  The  nature  of 
this  difference  must  be  considered,  and  the  DFT  results  must 


be  corrected. 


A. 2.  Deviation  of  the  DFT  from  the  Fourier  Transfrom 

In  appendix  C  it  is  shown  that  the  DFT  of  a  piecewise 
constant  field  deviates  from  its  continuous  transform  as 

DFm(x,y)]uu  =  kxuA _ kwA  (4>3)f 

^kxM'kyu)  4  sin[kx//2]  sin[kyuA/2] 

The  field  +(x,y,z),  is  assumed  to  be  adequately  represented 

by  a  piecewise  constant  field  with  the  field  value  of  each 

piece  given  by  one  corresponding  hologram  datum  value. 

Equation  (4.3)  is  therefore  relevant,  and  it  indicates  that 

while  the  DFT  of  the  hologram  data  will  provide  reasonable 

estimates  of  the  transform  field  for  small  k„  and  k„  ,  it 

y v 

will  overestimate  the  higher  spatial  frequency  region. 

Compounding  any  error  introduced  by  the  DFT  or  noise  in 
the  input  data  is  the  exponential  growth  of  the  reconstruc¬ 
tion  kernel  G”1  (z)  in  equation  (2.27)  with  spatial 

frequency  and  reconstruction  distance.  This  problem 
suggests  two  things:  First,  some  form  of  filtering  in 
k-space  will  be  necessary;  and  second,  deconvolutions  or 
source  field  reconstructions  which  include  the  evanescent 
wave  information  must  be  restricted  to  small  z  distances. 
Undersamping  in  k-space  is  another  source  of  error 

which  may  manifest  itself  if  either  4»(k  ,kv,z  )  or 


GI1^(kv'k,r>z)  varies  rapidly  over  the  DFT  sampling  intervals 
Ak^x  and  Aky^.  Undersampling  can  be  avoided  by  reducing  the 

k-space  sampling  intervals.  When  using  the  FFT,  this  sample 
spacing  is  most  easily  reduced  by  augmenting  the  hologram 
data  f(p,q,z  )  with  zeroes  in  the  same  manner  as  in 

t  (/, m)  of  equation  (3.7)  and  performing  the  FFT  of  the 

2N  x  2N  zero  augmented  hologram  array.  This  extended  array 
is  labeled  t'(p,q,z  )  in  analogy  with  t'  (i, m) . 

A. 3.  A  smoothing  function 


Based  on  the  above  considerations,  equation  (2.27)  can 
be  approximated  as  follows: 


«fD,N(/.m)=IDFT(W(kx)1,kyu)-DFT(H"(p,q,ZH)]-G^N(kx>l,kyv,ZH)) 

(4.4) 


where  W(kx^,kyy)  is  a  suitable  filter  funtion  and  t'(p,q,zH) 

is  the  hologram  data  augmented  with  zeroes.  For  the 
reconstructions  presented  later  in  this  paper,  the  following 
filter  function  was  used: 


W(kxji,kyv^ 


g(kr/kc-  l)/o 
z 

e(l  -  kr/kc)/a 


k_  <  k_ 
r  -  c 


(4.5) 


>  k_ 
r  c 


where  k*  =  k£  +  k^  and  k  and  a  are  adjustable  parameters 


of  the  filter.  Experience  has  shown  that  a  value  of  kc 

equal  to  0.6  (Nn/L)  and  a  value  of  a  equal  to  0.2  yield 
consistently  good  results.  In  a  particular  experiment  it 
may  be  desirable  to  change  these  parameters,  or  possibly  the 
filter  function's  form,  based  on  some  a  priori  knowledge  of 
the  source.  However,  any  filtered  source  field  reconstruc¬ 
tions  presented  or  compared  in  this  work  have  been  computed 
with  equation  (4.5)  and  the  filter  parameters  kc=0 . 6 ( Nn/L ) 

and  a  =0.2. 

B.  Reconstruction  and  inverse  transformation  test  results 


The  suitability  of  the  reciprocals  of  all  six  Green's 

function  forms  as  representatives  of  G-1(k  ,k  ,z)  in 

x  y 

equation  (4.4)  was  tested.  The  results  are  shown  in  the 
charts  of  figures  4.1a  and  4.1b  for  the  case  of  determining 

the  pressure  fields  P ( z0=0 • °65* )  and  z-component  of 
velocity  fields  v^  (  zQ=0 . 065A )  using  p°m(z=0.13A)  as  the 
input  field.  The  integrated  k-space  form  g^6^(z)  is  not 

0  t  N 


included  in  these  charts  as  the  errors  produced  using  this 
form  are  so  large  that  scaling  the  figures  to  accommodate 
them  obscures  any  information  about  the  other  forms.  The 
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Fig.  4.1  The  normalized  error  encountered  with  each  method 
in  the  source  field  reconstruction  of  (a)  the  pressure 
field  0.065  X  from  the  source  from  an  input  sample  array  of 
the  theoretical  pressure  field  0.13  X  from  the  source  and 
(b)  the  z-component  of  the  particle  velocity  field  0.065  X 
from  the  source  from  the  pressure  field  0.13  X  from  the 
source . 


sampled  k-space  Green's  function  g:  i(z)  is  the  best 

D  i  N 

representative  of  G"1  (k_,,k„,z)  in  all  reconstruction  cases 
considered . 

As  discussed  in  section  A. ,  some  form  of  filtering  must 
be  included  in  the  reconstruction  process  using  DFT ' s .  This 
need  for  filtering,  the  form  of  the  filter  function,  and  its 
effect  on  a  numerical  reconstruction  are  illustrated  in 
figure  4.2.  The  error  from  the  high  spatial  frequency 
approximation  associated  with  the  DFT  and  its  amplification 
by  the  reconstruction  process  is  apparent  if  the  DFT  of  an 

unfiltered  solution,  p^(z),  is  plotted  with  the  DFT  of 
p“m(z).  In  figure  4.2a  the  crosses  indicate  the  amplitude 
of  an  unfiltered  DFT  of  the  pressure  field  p^^ ( zo=0 . 065X ) 
as  processed  from  the  input  pressure  field  p°m( z=0 . 13X) . 

The  solid  line  plots  the  DFT  of  P/m ( z0=0 . 065X ) .  The 

difference  between  the  two  plots  clearly  grows  for  larger 
magnitudes  of  spatial  frequency.  In  real-space,  this  high 
frequency  error  shows  up  as  large  oscillations  of  the 
numerical  solution  about  the  actual  solution.  This  behavior 
is  evident  in  in  figure  4.2b  which  compares  the  unfiltered 

P/m^ ( zo=0 ’ 065X)  with  P;m( zo=0 . 065X ) .  The  numerical  solution 
is  vastly  improved  by  a  multiplication  with  the  filter 
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Fig.  4.2  (a)  plus  signs,  amplitude  of  the  DFT  of  a  sample  of  the  theoretical 

pressure  field  0.13  A  from  the  source  multiplied  by  the  [g^  Mz)]'1  to 

reconstruct  the  pressure  at  the  source.  The  solid  line  plots  the  theoretical 
values.  (b)  Theunf  ilte-red  field  against  theory,  (c)  the  filter  function  for 
three  values  of  a,  (d)  the  filtered  output  using  the  a  =  0 . 2  filter. 


function,  equation  (4.5).  This  function  is  plotted  in 
figure  4.2c  for  kc=  0.6  Nn/L  and  for  three  different  values 

of  a.  Values  of  kc=  0.6  Nn/L  and  a  =  0.2  have  been  found 

empirically  to  give  good  results.  These  values  are  used  in 
all  filtered  reconstructions  in  this  study  including  those 
which  provide  data  for  figures  4.1a  and  4.1b.  Figure  4 . 2d 
shows  the  excellent  fit  of  the  filtered  and  actual 
solutions . 

As  further  illustration  of  the  capabilities  of  the 
sampled  k-space  Green's  function  with  filtering,  figure  4.3a 

shows  a  comparison  plot  of  the  real  part  o  f  v^)(Zo=0.065X), 

where  the  input  field  is  p^m(z=0 . 13X) ,  with  v^m(zo=0 . 065X) 
plotted  as  the  solid  line;  figure  4.3b  shows  the  real  part 
of  v7m^zo=°)'  where  the  input  field  is  P°m( z=0 )  -  against 

v°m(zo=0).  Both  figures  show  good  agreement  between  the 

numerical  and  actual  solutions.  These  plots  along  with 
figures  4.1  and  4.2  indicate  that  accurate  solutions  to 
reconstruction  and  inverse  transformation  problems  are 
attainable  using  the  sampled  k-space  Green's  function, 

gi 5i(z)<  and  the  filter  function  of  equation  (4.5). 

For  reconstructions,  the  discussion  of  section  A.  and 
the  test  results  presented  in  this  section  indicate  that  the 

reciprocal  of  the  sampled  k-space  Green's  function  g^5^(z) 
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Fig.  4.3  (a)  Real  part  of  the  reconstructed  z-component  of 

the  particle  velocity  field  at  0.06  A  from  the  source 
starting  from  the  pressure  0.13  A  from  the  source.  (b)  Real 
part  of  the  z-component  of  the  particle  velocity  as 
reconstructed  in  the  source  plane  from  the  input  source 

plane  pressure.  The  reciprocal  of  g^  {z)  with  filtering 

was  used  in  both,  and  in  both,  the  solid  line  is  theory. 


with  a  spatial  frequency  filter  is  the  best  inverse 
propagator.  The  new  sampling  forms  are  of  little  use  with 
this  numerical  source  field  reconstruction  method. 


C.  Iterative  deconvolution 

For  practical  purposes,  it  is  not  likely  that  a  direct 
inversion  of  equation  (4.1)  would  ever  be  useful  due  to  the 
typical  size  of  the  matrix  involved,  4096x4096  in  this  work, 
and  due  to  its  nearly  singular  nature.  The  nearly  singular 
nature  of  this  matrix  follows  from  the  evanescent  decay  of 
high  spatial  frequency  input  fields  4*D  N  and  the  exponential 

enhancement  of  input  fields  with  spatial  frequencies  near 

the  propagating  wavenumber;  the  matrix  must  display  this 
decay  and  enhancement  if  it  is  to  model  wavefront 
propagation  accurately.  However,  the  fact  that  the  matrix 
multiplication  of  equation  (4.1)  can  be  carried  out  quickly 
with  the  FFT  suggests  the  search  for  an  iterative 
deconvolution  method  requiring  matrix  multiplications  and  no 
matrix  inversions. 

C.l.  Conjugate  gradient  descent  method 


reconstructions.  The  advantages  of  this  method  are  its 
objective  nature,  no  filter  parameters,  and  its  relative 
stability  in  comparison  with  the  filtered  solutions 

37 

presented  above.  Luenberger  describes  this  method  in 
terms  of  minimizing  the  functional 

f (x)  =  (£|q3)  -  2(t\x)  (4.6)  , 

where  Q  is  a  positive  definite  operator,  and  an  inner 

product  is  represented  by  (x|y).  The  solution  vector  xq  = 

Q-1B  uniquely  minimizes  this  functional.  Before  describing 
the  iterative  solution  of  this  problem,  the  correspondence 

is  indicated  between  x,  Q,  and  1$  and  the  acoustic  items 
+  .  the  Green's  function  matrix  ,(z  )]  of  equation 

(4.1),  and  the  hologram  field  t(z  ). 

H 

The  solution  vector  x  is  identified  directly  as  the 
collection  of  values  f/,x n)  representing  the  sought- 

after  boundary  condition.  The  positve  definite  operator  Q 
can  be  identified  as  the  Green's  function  matrix  premulti¬ 
plied  by  its  conjugate  transpose.  The  vector  B  can  then  be 
identified  as  the  result  of  premultiplying  the  collection  of 
values  t(p,q,zH)  by  the  conjugate  transpose  (indicated  by 

the  superscript  dagger, t)  of  the  Green's  function  matrix. 

The  positive  definite  quality  of  this  choice  for  the 
operator  Q  follows  from  the  observations  that  the  inner 
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product  (x|Q^)=x*Qx  is  then  the  norm  of  a  sample  of  values 
representing  a  radiating  wavefront  and  is  therefore 
positive,  and  the  only  bounded  field  which  can  be  zero 
over  a  plane  above  the  boundary  is  zero  everywhere  including 
the  boundary.  With  the  equation  (4.30) 's  components  thus 
identified,  it  is  clear  that  performing  a  source  field 
reconstruction  is  equivalent  to  minimizing  equation  (4.6). 

The  basic  idea  of  the  conjugate  gradient  descent  method 

is  to  move  from  each  successive  approximation,  given  by  xn, 

of  xq  in  a  direction  which  lies  near  to  the  negative 

gradient  of  equation  (4.6),  but  which  is  also  perpendicular 
to  all  previous  changes.  In  terms  of  these  successive 

approxima-  tions  xn  and  mutually  perpendicluar  directions  of 
change  j$n,  equation  (4.6)  can  be  written  as 

f(*n>  =  f<*n-i+*n3n>  (4>7)' 

The  ?n's  are  readily  chosen  to  minimize  the  functional  at 

each  step,  and  it  can  be  shown  that  the  direction  vectors 
]&ncan  be  calculated  recurrsively . 37  An  efficient  algorithm 

exists  where  xn  in  theory  converges  to  xq  in  a  finite  number 
of  steps,  and  the  algorithm  requires  only  one  evaluation  of 

O  Q 

the  operator  at  each  step.  The  unique  aspect  of  the 
application  of  this  technique  to  source  field  reconstruc- 
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tions,  or  deconvolution,  is  that  the  N2xN2  matrix  Q  needed 
never  be  formed.  All  matrix  multiplications  by  Q  can  be 
carried  out  exactly  by  the  discrete  convolution  algorithm 


presented  in  chapter  III  with  the  Green's  functions  gj.^1  or 

D  >  N 


C.2.  Results  from  the  implementation  of  this  method 


Although  the  conjugate  gradient  method  converges  in 
theory  to  the  unique  solution  ,  it  will  not  do  so 

numerically  due  to  roundoff  error.  However,  it  is 
reasonable  to  assume  that  only  the  highest  non-radiating 
spatial  frequency  modes  are  missing  from  the  non-unique 
numerical  solution  after  a  finite  number  of  steps,  and, 
therefore,  it  is  not  necessary  to  obtain  the  exact  unique 
solution  for  most  purposes.  A  practical  solution  is 
obtained  by  stopping  the  iterative  process  whenever  the 

magnitude  of  the  remainder  |  Q5<n— |  falls  below  a  chosen 
level . 

The  hidden-line  plot  in  figure  4.4  is  the  smoothed 
reconstruction,  by  the  method  of  section  IV. A.,  of  the  real 
part  of  the  normal  component  of  velocity  at  the  baffled 

piston  source,  v l ^ ( z  =0),  from  the  pressure  field  0.13X 


away,  p"m( z=0 . 13X ) .  Figure  4.5  shows  the  reconstruction 

obtained  after  43  iterations  by  the  techniques  of  this 

section  starting  with  the  same  input  field  p'jm(z=0 . 13X)  . 

The  standard  deviation  of  the  smoothed  reconstruction  from 
the  actual  boundary  field  is  more  than  five  times  that  of 
the  iterative  solution,  and  the  average  difference  is  nine 
times  that  of  the  iterative  solution. 

Figure  4.6  shows  the  solution  after  192  iterations. 

The  standard  deviation  of  this  solution  from  the  actual 
field  has  increased  by  15%  over  that  obtained  after  43 
iterations  inspite  of  the  fact  that  the  magnitude  of  the 

remainder  |Qxn-£|  has  decreased  by  a  factor  of  10  from 

0.00l|$|  to  0 . 000 1 | $ | .  The  important  point  here  is  that 
the  first  40  iterations  reduced  the  remainder  by  a  factor 
of  1000,  while  nearly  five  times  as  many  iterations  yielded 
only  a  tenfold  reduction  and  did  not  improve  the  solution. 
This  suggests  that  the  iteration  process  may  be  objectively 
stopped  when  it  is  observed  that  the  magnitude  of  the 
remainder  is  no  longer  decreasing  rapdily. 

The  major  disadvantage  to  this  technique  is  the  amount 
of  time  required;  nearly  ten  times  as  much  total  job  time  is 
required  to  compute  a  solution  by  40  iterations  as  is 
required  for  a  smoothed  solution.  The  independence  from 


Fig.  4.5  The  reconstruction  obtained  of  the  real  part  of 
the  normal  velocity  component  at  the  baffled  piston  source 

v/m^(z0=0)»  from  the  pressure  field  0.13X  away  after  43 

iterations  by  the  techniques  of  section  IV. C. 


adjustable  filter  factors  and  the  possibility  of  increased 
accuracy  may  warrant  such  an  increase  in  computation  time. 
Further,  the  actual  computation  time  may  be  reduced  by 
performing  the  calculations  on  an  array  processor  and  by 
improving  the  algorithm  itself.  The  execution  time  of  the 
iterative  solution  is  likely  to  improve  much  more  than  the 
smoothed  solution  with  the  use  of  an  array  oriented  machine 
since  those  sections  which  are  least  ammenable  to  array 
manipulation,  intialization  and  input/output,  are  virtually 
identical  in  both  cases  and  constitue  the  bulk  of  the 
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latter.  The  basic  algorithm  can  be  generalized  so  that 
each  step  requires  only  one  matrix  multiplication  by  [Gd  n] 

and  not  the  two  matrix  multiplications,  [G*  1  and  [ G  1, 

represented  by  Q;  this  would  cut  the  iteration  time  in  half. 
The  conjugate  gradient  descent  method  using  the  FFT  to  carry 

out  discrete  convolutions  with  the  Green's  functions  g^2^ 

and  gi4i,  as  described  in  chapter  III,  deserves  further 
study  in  actual  Nearfield  Acoustic  Holography  experiments. 


Chapter  V. 


SOLUTION  OF  THE  SURFACE  PROBLEM  FOR  ODD-SHAPED  RADIATORS 

In  many  practical  acoustic  problems  with  odd-shaped 
radiating  surfaces,  the  normal  derivative  of  the  field  in 
the  form  of  the  normal  component  of  the  surface  velocity  is 
specified,  and  the  pressure  field  over  the  surface  is 
sought.  Forming  a  numerical  solution  to  this  problem  is 
also  the  first  step  in  the  reconstruction  process  covered  in 
chapter  VI.  Sections  A  and  B  of  this  chapter  develop  a 
numerical  solution  for  this  problem  and  present  the  theory 
behind  this  method.  In  section  C,  the  non-trivial  but 
theoreticaly  tractable  spherical  source  boundary  is 
presented  as  a  suitable  test  case  for  the  techniques  of  this 
chapter  and  the  next . 


A.  Relating  the  surface  velocitv  and  pressure 


For  a  radiating  surface  which  does  not  correspond  to  a 
level  surface  in  a  separable  coordinate  system  of  the 
Helmholtz  equation,  there  is  no  simple  integral  which  yields 
the  field  4*  in  terms  of  either  the  Dirichlet  boundary 


condition  +  or  the  Neumann  boundary  condition  -19+  alone 

9n 


specified  over  the  surface  .  With  the  requirement  that  the 
field  satisfy  radiation  conditions,  the  specification  of 
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either  Dirichlet  or  Neumann  boundary  conditions  is 

O  Q 

sufficient  to  determine  +  uniquely  for  any  closed  surface^  . 
The  Surface  Helmholtz  Integral  Equation  and  the  Interior 
Helmholtz  Integral  Equation  relate  the  two  boundary 
conditions  accordingly.  The  use  of  these  equations  in 
continuous  form  and  in  discrete  approximation  are  discussed 
in  this  section  while  their  actual  reduction  to  discrete 
form  is  discussed  in  the  next  section.  The  last  section  of 
this  chapter  will  present  test  results  for  this  method. 

A . 1 .  The  Surface  Helmholtz  Integral  Equation 


In  chapter  II,  the  Surface  Helmholtz  Integral,  equation 
(2.19),  was  derived  as  a  limiting  case  of  the  Helmholtz 
Integral  for  the  general  field  +  .  Written  in  terms  of  the 

acoustic  pressure  field  p(rg)  and  the  surface  normal 


component  of  particle  velocity  vn(rg),  this  equation  is 


22  f (?e-rl)dS 


ipckO  vn ( rg )  Gf(rg-r^)dS 

J  g 


(5.1) . 

The  integrations  in  this  equation  must  be  replaced  by 
appropriate  numerical  integration  approximations  to  solve 
the  general  problem  of  odd-shaped  radiators.  At  some  point 
in  their  application,  most  numerical  integration  schemes 


reduce  to  matrix  operations  on  a  set  of  coefficients 
representing  the  known  and  unknown  fields  vn(rs)  and  p(rs) . 

Denoting  these  coefficient  sets  as  ^(^g)^  and  pCr)^  the 

matrix  approximation  to  equation  (5.1)  is 

£i|s>q  =  DSqj  P<rs)j  "  MSqj  VVj  (5-2)' 

with  the  Einstein  summation  convention  in  force.  Equation 
(5.2)  can  be  written  as 

C5q;|  p(?s)j  =  MSqj  vn(?)j  (5.3) 

where  =  ^qj  ~  -j*qj  anc*  't^ie  matrices  [ D ]  and  [M]  are 

appropriate  to  the  chosen  numerical  integration  scheme. 

At  certain  characteristic  frequencies,  equation  (5.1) 
is  not  sufficient  to  determine  a  unique  solution  to  the 

radiation  problem . 30 • 39 ' 40  The  problem  is  that,  at  these 
frequencies,  interior  eigenmodes  exist  which  are  homogenous 
solutions  over  the  surface  under  study  •  To  any  particular 
solution  of  the  problem,  it  is  then  possible  to  add  one  or 
more  homogeneous  solution.  This  problem  manifests  itself  as 

nearly  singular  matrices  [E5]  and  [MS]  in  the  approximation 
of  equation  (5.3)  . 


The  Interior  Helmholtz  Integral  Equation 


The  homogeneous  solutions  of  the  Helmholtz  equation  are 
not  radiating  solutions,  and  it  is  necessary  to  incorporate 
the  Interior  Helmholtz  Integral  Equation,  equation  (2.18), 

on 

with  the  surface  equation  to  eliminate  these  solutions'^ . 
Written  in  terms  of  the  acoustic  pressure  field  p(rs)  and 


the  surface  normal  component  of  particle  velocity  vn(rs)< 

the  Interior  Helmholtz  Integral  Equation  is 
r>  p 


0=<p  p( )  ggffrj-r^JdS  -  ipckfj)  vn(r^)  Gf(rI-r,l)dS  (5.4) 


where  the  point  is  interior  to  the  surface.  Equation 


(5.4)  can  be  approximated  by 

DIqj  =MIqj  V?»j  <5-5> 

where  the  subscript  q's  refer  to  a  number  of  points  rj  in 

the  interior.  Equations  (5.3)  and  (5  5)  can  be  treated  as  an 
overdetermined  system  of  equations  and  solved  accordingly. 

In  the  literature , ^ ' 40  a  largely  unanswered  question  is  how 
to  choose  the  interior  points,  as  points  which  have  nearly 
zero  field  with  an  interior  eigenmode  boundary  solution  must 


be  avoided.  The  next  section  details  the  reduction  of 
equation  (5.1)  to  equation  (5.2)  and  also  introduces  a 


method  for  detecting  characteristic  frequencies  and  for 
choosing  interior  points  to  avoid  the  problems  caused  by 
these  solutions. 


B.  Reduction  to  finite  and  discrete  operations 

The  reduction  of  the  surface  problem  to  the  matrix  form 
of  equations  (5.3)  and  (5.5)  and  their  solution  is  presented 
in  three  parts  in  this  section.  The  geometric  represen¬ 
tation  of  the  surface  as  an  assembly  of  flat  triangular 
plates  along  with  the  corresponding  scheme  for  performing 
the  numerical  integrations  is  presented  first.  The 
detection  of  characteristic  frequencies  for  odd-shaped 
radiators  is  covered  next.  The  formation  of  a  surface 
solution  matrix  is  considered  last. 

B.l.  Finite  elements 


The  continuous  surface  of  an  odd-shaped  radiator  must 
be  represented  by  a  generally  simpler  geometric  surface 
which  can  be  specified  by  a  finite  number  of  parameters. 

The  method  of  representing  the  surface  presented  here  is 
essentially  the  method  most  recently  described  by  Koopmann 

and  Benner.40  The  actual  surface  is  replaced  by  a  number 
of  interconnected  flat  triangular  plates  whose  vertices  lie 
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on  the  actual  surface.  The  acoustic  pressure  and  normal 
component  of  particle  velocity  are  considered  constant  over 
each  triangular  element  and  equal  to  their  respective  values 
at  the  center  of  each  element.  Since  the  elements  are 
flat,  their  normal  directions  are  easily  determined  by 
transforming  the  global  set  of  coordinates  in  which  the 
problem  is  originally  specified  into  new  sets  of  coordinates 
for  each  element.  These  local  sets  of  coordinates  are 
defined  by  the  three  element  vertices  in  such  a  way  that  the 
local  z-axis  points  in  the  local  normal  direction. 

By  approximating  the  surface  and  surface  field  with 
these  triangular  elements,  equation  (5.1)  can  be 
approximated  as 

1 

£(?_)=  Z  P ( r j )  9Gf (r  -r' )dx'dy' 

j*q  3  3zr  q 
j 


-  ipck  Z  vn(rj)  Gf (rq-r ' )dx'dy'  (5 


where  x,  y,  and  z  are  the  local  coordinates  for  each 


trianglular  element  j  and  the  vectors  r^  locate  the  center 

of  each  triangle  globally.  The  integrands  are  simply  half 
the  Green's  functions  GD  and  GN  encountered  in  chapters  II, 

III,  and  IV.  The  first  integral  for  j=q  yields  zero 
because  the  element  is  flat  and  the  entire  integral  of 


equation  (5.1)  omits  a  small  region  around  the  point  r^; 


thus,  integration  of  GD  over  element  q  is  omitted  in 
equation  (5.6) . 

By  multiplying  equation  (5.6)  by  two,  substituting  GD 
and  Gn  for  the  integrands,  noting  that  integrating  GD  over 

the  entire  element  q  yields  one,  and  finally  rearranging  the 
integrals  with  respect  to  the  equals  sign,  equation  (5.6) 
can  be  written  as 


P(rj> 


GD(r  )dx'dy ' 


*  ipck  E  vn( r j ) 


G„(*  -?')dx'dy'  (5.7). 


The  integrals  are  the  same  type  as  in  equation  (3.5) 
differing  only  in  that  the  integration  areas  here  are 
triangles  as  opposed  to  squares.  The  integrations  may  be 
approximated  by  expanding  the  Green's  functions  as  is  done 
in  chapter  III.  The  details  of  this  expansion  and  the 
resulting  quadrature  formulae  are  presented  in  Appendix  B. 
The  matrix  coefficients  of  equation  (5.2)  can  then  be 


identified  as 


»■  V  T  V  V  l.*  V  »  V  ',  ’■V  -.'  »■■  rj  rr.yjn'»-;  vu  r: 


my* V9imy\*y.»i'y>wrv*jvm'myim  1 1 
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points  if  solutions  at  these  frequencies  are  needed.  In  this 
section,  a  method  for  detecting  these  frequencies  is 
presented  and  a  suggestion  on  choosing  the  interior  points 
is  developed. 

It  is  assumed  that  the  surface  equations  can  be 
adequately  represented  by  a  set  of  linear  equations,  as  in 
equation  (5.3),  throughout  the  frequency  region  of  interest. 
Therefore,  if  the  continuous  equation  admits  non-radiating 
homogeneous  solutions,  its  numerical  representation  as  a 
linear  system  of  equations  should  admit  these  solutions  as 
well.  Exactly  how  these  non-radiating  interior  solutions 
find  their  way  into  an  exterior  approximation  formulation  is 
intuiatively  clear  if  the  interior  problem  for  the  same 
surface  is  considered. 

The  monopole  matrix  [MS]  of  equation  (5.3)  for  the 
interior  problem  differs  from  that  for  the  exterior  problem 
only  by  a  minus  sign  for  each  element,  while  the  dipole 


matrices  [E3]  are  identical.  At  characteristic  frequencies, 
velocity  eigenmodes  exists  for  the  interior  which  yield  zero 
surface  pressure  by  definition.  The  left-hand  side  of 
equation  (5.3)  is  zero  for  such  a  mode  and  so  reduces  to 

MSqj  vi?<?S5=  -MSq j  vn^s>=  °'  v£<*s>*  0  <5-10>' 

where  the  superscripts  I  and  E  indicate  the  matrices  of  the 
interior  and  exterior  problems  and  where  the  superscript  0 
identifies  the  vector  as  representative  of  an  eigenmode.  If 


equation  (5.10)  has  a  solution,  then  the  equation 

p°(rs)=  BSjjj  P°(rs)=  0,  p°(rs)*  0  (5.11) 

must  have  a  solution  p°(rs),  which  is  the  radiating  solution 
corresponding  to  -vn(rg).  Equations  (5.10)  and  (5.11)  can 

only  be  satisfied  if  the  matrices  [MS]  and  [175]  are 

singular41;  that  is, their  determinants  are  zero.  The 
problem  of  determining  the  characteristic  frequencies  of  an 
odd-shaped  surface  is  therefore  equivalent  to  measuring  the 
relative  singularity  of  the  matrices  as  a  function  of 
frequency. 

According  to  Stewart,41  the  difficulty  in  determining 
numerical  singularities  is  that  the  matrices  involved  are 
only  approximations  to  the  actual  problem  under  study  and  so 
may  only  be  close  to  singular  matrices  when  the  real  problem 
is  singular.  Furthermore,  many  methods  for  decomposing  a 
matrix  may  transform  a  nearly  singular  matrix  to  a  form 
which  is  not  nearly  singular.  However,  the  Singular  Value 
Decomposition  yields  a  nearly  singular  matrix  if  and  only  if 

the  original  matrix  is  nearly  singular.41  Therefore,  this 
decomposition  is  ideally  suited  for  determining  the 
characteristic  frequencies  of  an  odd-shaped  surface  and  the 
eigenmodes  as  well.  To  show  how  this  is  possible,  the 
Singular  Value  Decompostion  is  examined. 


By  premultiplying  and  postmultiplying  a  general  matrix 
by  a  set  of  unitary  transformations,  the  Singular  Value 
Decomposition  reduces  the  matrix  to  diagonal  form.  By 

gathering  these  transf romations  into  the  matrices  [U]  and 

[V*]  (the  dagger  indicates  the  conjugate  transpose),  the 
general  mxn  matrix  [X] ,  where  m>n,  can  be  written  as 

[X]  =  [U] [§]  [Vf]  (5.12)  , 

where  [U][U*]=[I],  [ V] [ V* ]=[ I ] ,  and  [E]  is  an  nxn  diagonal 
matrix  with  diagonal  elements  8^,  the  so-called  singular 
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values.  It  can  be  readily  shown  that  the  singular  values 
measure  the  relative  output  vector  magnitude  for  each  of  a 
complete  set  of  input  vectors,  namely  the  columns  of  (V]. 

(  For  a  hermitian  matrix,  the  singular  values  are  the 
eigenvalues,  and  the  columns  of  (V]  are  the  eigenvectors.) 

In  decomposing  either  [MS]  or  [DS] ,  this  set  covers  all 
possible  surface  velocity  or  pressure  distribution 
approximations.  It  is  the  low  frequency,  long  wavelength 
regime  which  is  of  main  interest  in  this  work,  and  it  is 
assumed  that,  if  interior  eigenmodes  exist,  they  must  be  few 
in  number  relative  to  the  number  of  degrees  of  freedom  of 
the  approximation  scheme.  Since  all  of  the  columns  of  [V] 
are  orthogonal,  eigenmodes,  if  they  exist,  must  be 
represented  by  a  distinct  set  of  columns  in  [V]  possessing 
correspondingly  small  singular  values  in  comparison  with  the 
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singular  values  of  the  majority  of  other  modes.  Therefore, 
the  presence  of  a  very  small  singular  value,  relative  to  the 
remainder  of  values  for  a  given  matrix,  indicates  the 
existence  of  an  interior  eigenmode,  at  least  in  a  numerical 
sense . 

The  lower  characteristic  frequencies  of  an  odd-shaped 
surface  enclosed  volume  can  be  found  by  forming  and 

decomposing  the  matrices  [MS]  and  [D5]  over  a  range  of 
frequencies.  In  a  plot  of  the  calculated  minimum  singular 
value  ratio  as  a  function  of  frequency,  these  frequencies 
should  be  characterized  by  a  pronounced  dip  in  the  plot. 

An  example  of  this  is  presented  in  the  last  section  of  this 
chapter.  Notice  that  only  the  singular  values  need  be 
calculated  for  this  process,  and  this  can  be  done  in  much 
less  time  then  is  required  for  the  complete  decomposition 
as  in  the  form  of  equation  (5.12). 

When  a  solution  is  required  at  an  identified  character¬ 
istic  frequency,  the  formation  of  the  matrix  [V]  of  equation 
(5.12)  is  useful.  Since  a  column  q  of  [V]  which  is 
associated  with  a  relatively  small  singular  value  S_  is  an 

'■d 

eigenmode  approximating  vector,  the  suitability  of  a  set  of 
interior  points,  to  be  used  in  forming  equation  (5.5),  can 
then  be  expected  if  the  following  conditions  are  met: 

DI  p°(?_)_  >>  0  and  MImri  v°  (?_)_  >>  0  (5.13), 

mn  s  n  mn  nq '  s '  n  '  '  ' 
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where  pG  and  vG  are  all  the  q^*1  columns  of  the  [V]  matrices 
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obtained  in  decomposing  respectively  [E15]  and  [MS]  such  that 

Sg  is  very  small.  Equation  (5.13)  insures  that  the  set  of 

chosen  interior  points,  along  with  the  Interior  Helmholtz 
Integral  Equation,  yield  a  system  of  equations  which  are  not 
satisfied  by  the  non-radiating  boundary  conditions. 


B.3.  Formation  of  the  surface  solution  matrix 


The  goal  of  this  chapter  is  to  form  an  invertible 
matrix  which  relates  the  surface  pressure  and  normal 
component  of  surface  velocity  as  follows: 

P<*S>q  =  GSqj  V?S>j  (5.14). 

At  non-characteristic  frequencies,  the  matrices  of  equation 
(5.3)  should  be  non-singular,  and  the  matrix  [GS]  may  be 
formed  directly  as 


GSqj  =  MSnj 

At  characteristic  frequencies,  the  matrix  must  be 
generated  from  the  overdetermined  system  of  equation 
and  equation  (5.5).  The  overdetermined  system  may  be 
written  in  the  form  of  a  single  matrix  equation  as 


(5.15). 


(5.3) 


U5 

-DI 


P(?s) 


MS 

MI 


V?s> 


(5.16)  . 


are  unitary  matrices.  The  matrix  [  GS  ]  is  then  given  by  the 
expression 


[  GS  ]  =  [  V  ][  S_1][  Uj4 


W. 


1 1 
?21 


[  r  ][  zf  ] 


(5.21)  . 


C.  Testing  the  techninque 


To  test  the  suitability  of  the  techniques  described  in 
this  chapter,  and  those  to  be  described  in  the  next  chapter, 
a  closed  boundary  for  which  non-trivial  but  theoretically 
tractable  fields  exists  is  desirable.  As  in  so  many 
problems  in  physics,  the  sphere  is  an  ideal  test  case. 
Although  the  spherical  coordinate  system  is  one  of  the 
separable  coordinate  systems  of  the  Helmholtz  Equation, 
no  advantage  is  taken  of  this  symmetry  in  the  techniques  of 
this  chapter  or  the  next;  thus,  in  this  sense,  the  surface 
of  a  sphere  may  be  taken  as  an  odd-shaped  surface.  In  this 
section  the  approximation  of  the  surface  is  presented,  and 
the  characteristic  frequencies  predicted  with  this 
approximation  by  the  technique  of  section  V.B.2.  are 
compared  with  those  of  a  true  spherical  surface. 
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C.l.  Approximating  the  spherical  surface 

The  surface  of  the  sphere  is  approximated  in  this  work 
by  a  60-faced  polyhedron.  Each  of  the  faces  are  triangular, 
and  the  entire  set  is  described  by  32  vertex  or  nodal  points 
which  lie  on  a  circumscribed  sphere.  These  32  nodal  points 
are  the  20  vertices,  of  a  dodecahedron  plus  12  additional 
points  obtained  by  radially  projecting  the  center  of  each 
face  of  the  dodecahedron  onto  the  circumscribing  sphere. 
Field  values  are  calculated  at  the  center  point  of  each 
triangular  face.  These  sixty  points  all  lie  on  an  inscribed 
sphere  whose  radius  is  0.923  times  that  of  the 
circumscribing  sphere.  This  discrepancy  must  be  accounted 
for  in  any  comparison  with  theory  for  the  sphere.  The  author 
was  spared  the  task  of  calculating  the  cartesian  coordinates 
of  these  nodal  points  and  of  determining  which  points 
describe  each  triangle  through  the  efforts  of  an  industrious 
undergraduate,  Ray  Sova,  who  researched  this  particular 

A  5 

problem.  The  polyhedron  is  pictured  in  figure  5.1.  The 
edges  of  each  adjacent  triangle  are  separated  for  graphical 
purposes  only. 
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C.2.  A  comparison  of  theoretical  and  experimental  values 
for  the  characteristic  frequencies  of  a  sphere 

The  monopole  and  dipole  surface  matrices  of  equation 
(5.3)  where  calculated  for  the  polyhedral  representation  of 
the  sphere  over  a  range  of  wavenumbers  from  near  zero  up  to 
twelve.  In  all  cases  the  radius  of  the  circumscribing  sphere 
was  fixed  at  unity.  The  minimum  singular  value  ratio  of  the 
monopole  matrix  as  a  function  of  wavenumber  is  plotted  in 
figure  5.2. 

The  spherical  Bessel  functions  are  the  solution  to  the 
radial  part  of  the  Helmholtz  Equation  in  spherical 

coordinates  for  the  interior  problem.46  Therefore, 
eigenmodes  or  resonances  occur  whenever  the  frequency 
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normalized  radius,  ka,  equals  a  zero  of  a  spherical  bessel 

function.  The  theoretical  resonances  or  characteristic 
frequencies  for  a  sphere  with  radius  equal  to  that  of  the 
inscribed  sphere  of  the  polyhedron  are  indicated  by  the  plot 
border's  outer  tick  marks. 

The  same  type  of  plot  is  presented  in  figure  (5.3)  for 
the  minimum  singular  value  ratio  of  the  dipole  matrix  of 
equation  (5.3).  Both  plots  show  excellent  agreement  between 
theory  and  calculation.  As  the  frequency  or  wavenumber 
increases,  this  agreement  should  break  down  as  more  complex 


eigenmodes  are  possible.  A  conservative  estimate  of  the 
useful  range  of  the  60-sided  polyhedron  for  representing  the 
sphere  is  obtained  by  equating  the  shortest  distance  between 
two  triangular  element's  centers  with  one  half  the  minimum 
valid  wavelength.  This  cutoff  falls  around  a  wavenumber  of 
five . 

To  produce  the  data  for  the  plots  of  figures  (5.4)  and 
(5.5),  a  number  of  interior  points  lere  used  to  eliminate 
the  interior  eigenmodes  as  possible  solutions  to  the 
numerical  formulation  of  .he  exterior  problem.  The  singular 
values  are  those  of  the  overdetermined  system  of  equation 
(5.18).  This  correction  was  performed  over  the  range  between 
the  dashed  lines.  The  elimination  of  the  numerical 
singularity  is  clear  in  these  plots. 
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Fig.  5.4  Monopole  matrix  characteristic  frequency  curve 
for  the  augmented  polyhedron  model.  (Augmented  between 
dashed  lines.) 
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Fig.  5.5  Dipole  matrix  characteristic  frequency  curve 
for  the  augmented  polyhedron  model.  (Augmented  between 
dashed  lines.) 


Chapter  VI . 


RECONSTRUCTION  OF  ODD-SHAPED  SOURCES 


The  extension  of  the  techniques  of  Nearfield  Acoustic 
Holography  to  include  the  reconstruction  of  odd-shaped 
sources  is  highly  desireable  for  both  theoretical  and 
practical  applications.  It  is  the  aim  of  this  chapter  to 
propose  a  numerical  method  for  accomplishing  this  extension 
and  to  present  the  results  of  simulation  experiments  which 
support  this  proposal . 

Generalized  Holography  in  any  separable  coordinate 
system  of  the  Helmholtz  Equation  is  based  on  the  analytic 
transformation  of  fields  to  and  from  an  appropriate  eigen- 
space  wherein  wave  propagation  is  performed  by  a  simple 

Q 

multiplication  operator.  As  is  shown  earlier,  in  the  case 
of  planar  holography,  source  field  reconstruction  is 
accomplished  by  Fourier  transforming  the  measured  field  and 
then  dividing  the  transformed  field  by  the  appropriate 
Green's  function  of  equation  (2.25  or  2.26).  In  chapter  IV 
it  was  shown  how  this  type  of  transformation  approach  is 
highly  desireable  In  reconstruction  problems  as  experimental 
or  computational  noise  tends  to  be  transformed  to  the  high 
order  modes  of  the  eigenspace  and  not  to  the  lower  modes. 


Being  an  inversion  process,  source  field  reconstruction 
amplifies  any  data  noise.  By  relagating  the  noise  to  higher 
order  modes  with  an  appropriate  eigenspace  transformation,  a 
reconstruction  of  the  lower  order  modes  is  possible. 

An  analytic  transformation  is  not  possible  in  the  more 
general  case  of  odd-shaped  source  boundaries.  However,  it  is 
possible  to  generate  and  perform  the  appropriate  trans¬ 
formations  numerically.  The  Singular  Value  Decomposition 
plays  a  central  role  in  this  general  reconstruction  process 
much  as  the  Fast  Fourier  Transform  does  in  the  case  of 
planar  holography.  Since  the  Singular  Value  Decompostion 
operates  on  a  matrix,  its  application  requires  the  explicit 
development  of  the  linear  relationship  between  the  normal 
component  of  surface  velocity  and  the  radiated  field 
pressure  and  a  reduction  of  this  relationship  to  a  finite 
matrix  form. 

A.  Relating  the  surface  velocity  to  the  field  pressure 


For  an  odd-shaped  radiator,  reconstruction  is  defined 
as  the  determination  of  the  normal  component  of  the  surface 
velocity  at  the  source  from  a  measurement  of  the  radiated 
pressure  field  at  a  number  of  points  above  its  surface. 
Experimentally,  these  measurements  may  be  obtained  with  a 
planar  array  of  microphones,  or  they  may  be  obtained  with  a 


more  flexible  microphone  system  at  a  number  of  points 
surrounding  the  surface.  A  general  method  for  forming  a 
matrix  relating  the  surface  velocity  to  the  radiated 
pressure  measured  at  these  microphones  is  a  fundamental 
requirement  of  the  reconstruction  process.  To  lay  the 
groundwork  for  such  a  method,  the  continuous  integral 
relationships  between  the  normal  component  of  surface 
velocity  and  the  radiated  pressure  in  the  field  are  examined 
in  this  section. 

A . 1  The  Helmholz  Integral 


The  Helmholz  Integral,  equation  (2.17),  can  be  written 
in  terms  of  the  acoustic  pressure  field  and  the  particle 
velocity  field  as: 


p(rf)=4>  p(  Tg )  9Gf  ( rf-rJ.)dS  -  ipcko  vn(r^)  Gf(rf-r^)dS 
J  e  dn  J  q 

a  b  (6.1) 


This  expression  relates  the  acoustic  pressure  at  a  point  r^ 

in  the  field  above  the  surface  S  to  the  surface  pressure  and 
the  sought-after  normal  component  of  surface  velocity.  The 
Surface  Helmholz  Integral  Equation  (5.1)  and  the  Interior 
Helmholtz  Integral  Equation  (5.4)  provide  the  link  between 
the  surface  pressure  and  velocity  as  discussed  in  the 
previous  chapter.  Together,  equations  (5.1),  (5.4),  and 


(6.1)  uniquely  determine  the  pressure  in  the  field  in  terms 
of  the  normal  component  of  surface  velocity. 


A. 2  Experimental  uniqueness  problems 

In  the  last  chapter,  a  fundemental  uniqueness  problem 
was  encountered  with  the  use  of  the  Surface  Helmholtz 
Integral  Equation  due  to  the  presence  of  homogeneous 
solutions  of  an  interior  problem  closely  related  to  the 
exterior  radiation  problem  under  study.  It  is  assumed  that 
the  techniques  of  the  last  chapter  are  sufficient  to 
overcome  this  problem.  With  the  use  of  equation  (6.1)  in 
this  chapter  however,  a  different  source  of  non-uniqueness 
is  encountered;  this  non-uniqueness  can  be  viewed  as  more 
experimental  in  nature  than  the  unavoidable  problem  with 
charateristic  frequencies  encountered  in  the  last  chapter. 

Given  the  normal  component  of  surface  velocity  for  the 
odd-shaped  source  and  the  Surface  and  Interior  Helmholtz 
Integral  Equations,  a  unique  surface  pressure  field  is 
determined.  These  two  surface  fields  in  turn  define  a 
unique  pressure  field  above  the  surface  by  way  of  equation 
(6.1).  However,  this  radiated  pressure  field  can  only  be 
observed  experimentally  over  a  limited  region  and  with 
finite  resolution.  In  analogy  with  the  theory  of  planar 
holography,  it  is  intuitive  that  spatial  frequency  modes  of 
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an  odd-shaped  surface  must  generate  evanescent  fields  when 
their  wavenumber  exceeds  that  of  the  propagating  wavenumber. 

The  information  about  these  modes  is  diminished  at  the  field 
measurement  loci  relative  to  the  information  from  lower 
frequency  radiating  modes.  If  the  evanescent  field  of  one 
mode  falls  below  the  dynamic  range  of  the  measuring  system 
relative  to  a  second  mode,  surface  modes  with  or  without 
this  evanescent  mode  will  appear  experimentally  identical. 

Numerically  this  problem  manifests  itself  as  a  singular 
system  of  linear  equations.  A  reconstruction  solution  to 
this  problem  is  possible,  but  it  requires  special  treatment 
of  all  information  about  these  modes. 

B.  Finite  and  discrete  operations 

% 

I 

The  numerical  reconstruction  problem  can  be  broken  down 
into  two  parts.  First,  a  single  matrix  is  formed  relating 

the  finite  set  of  numbers  vn(rs)j,  which  represent  the  | 

normal  component  of  surface  velocity,  to  a  measured  set  of 
field  pressure  values  p(rj)g.  And  second,  a  minimal  solution  ' 

to  this  matrix  problem  is  determined.  These  two  problems  are  ! 


considered  in  this  section. 


2 


B . 1  Formation  of  the  field  matrix 

In  order  to  reduce  the  size  of  the  numerical 
odd-shaped  surface  reconstruction  problem,  it  is  desireable 
to  form  a  single  matrix  which  directly  relates  the  normal 
component  of  surface  velocity  to  the  pressure  at  a  number 
of  points  in  the  field  above  the  surface.  This  matrix  is 
refered  to  as  the  field  matrix  [GF],  and  with  it  the 
following  equation  can  be  written: 

P<rf)q  =  GFqj  V*s>j  <6-2>- 

As  already  discussed  above,  there  is  no  single  integral 
relating  the  surface  velocity  and  field  pressure,  although 
the  pressure  is  uniquely  determined  by  the  system  of 
equations  (5.1),  (5.4),  and  (6.1).  The  reduction  of  equation 

(5.1)  and  (5.4)  to  finite  form  was  covered  in  the  last 
chapter,  and  the  same  methods  used  there  can  be  applied  to 
the  reduction  of  equation  (6.1)  to  a  finite  and  discrete 
form.  By  treating  the  odd-shaped  surface  as  constructed  of 
a  number  of  edge-connected  flat  triangular  elements,  each 
with  constant  surface  pressure  and  velocity,  as  in  chapter 
V,  equation  (6.1)  can  be  approximated  by  the  following 
equation : 

p(?f)q  =  DFqj  p(?s)3  +  MFgj  vn(?s)j  (6.3), 

where  the  matrix  coefficients  can  be  identified  as 


\  j  GD(xq-x' ,yq-y ' ,zq-z' )dx'dy'  (6.4) 


and  MF, 


=  ^  gn ( xq-x 1 , yq-y ' , zq-z ' ) dx '  i 

j  *  *  i 


6.5). 


These  integrals  differ  by  the  simple  constant  factors  of 
plus  and  minus  one-half  from  the  integrals  of  equations 
(5.8)  and  (5.9)  respectively.  Their  numerical  evaluation  is 
presented  in  appendix  B. 

By  using  the  surface  solution  matrix  [GS]  of  equation 
(5.21),  the  surface  pressure  vector  pC?s)j  can  be  written  in 
terms  of  the  surface  velocity.  With  this  replacement, 


equation  (6.3)  can  be  written  as 

P<?f>q  =  <DPqn  GSnj  +  MFq;J  >  V*s  >  j 
and  the  elements  of  [GF]  can  be  identified  as 


(6.6)  , 


GF  •  =  (DF  GS  •  +  MF  • ) 

qj  qn  ^Jnj  ljrqj  ’ 


(6.7)  . 


Equation  (6.7)  provides  the  formula  for  constructing  the 
elements  of  [GF]  from  a  set  of  readily  calculated 
coefficients  for  any  odd-shaped  geometry. 


B.2  Solution  with  the  Singular  Value  Decomposition 


The  generalized  inverse  of  a  matrix  is  defined  in 
equation  (5.18)  of  chapter  V  in  terms  of  the  matrix's 


Singular  Value  Decomposition.  The  generalized  concept  of  an 
inverse  is  needed  in  chapter  V  inorder  to  obtain  a  single 
nxn  matrix  relating  the  surface  velocity  and  surface 
pressure  from  the  mxn,  m>n,  overdetermined  system  of 
equations  (5.3)  and  (5.5),  where  n  is  the  number  of  degrees 
of  freedom  for  the  surface  model,  and  m-n  is  the  number  of 
interior  points.  In  this  chapter,  the  Singular  Value 
Decomposition  is  used  in  a  similar  way,  but  it  is  not 
necessary  to  actually  form  the  generalized  inverse  of  the 
field  matrix  [GF]  as  is  necessary  for  [GS]  in  chapter  V.  It 
is  interesting  and  instructive  to  clarify  these  ideas  by 
comparing  the  numerical  process  of  wavefront  propagation  and 
source  field  reconstruction  for  an  odd-shaped  source  with 
the  analytic  process  for  a  source  with  a  level  surface  in  a 
separable  coordinate  system. 

For  the  odd-shaped  source  problem,  the  surface 
is  represented  by  a  model  with  n  degrees  of  freedom  and 
the  field  above  the  surface  is  sampled  at  m  points  where  m 
is  generally  larger  than  n.  The  mxn  matrix  [GF]  then  has 
a  decomposition  with  the  same  form  as  equation  (5.17),  that 
is 
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(6.8)  . 


The  matrix  [GF]  is  the  operator  which  produces  the  pressure 
at  a  number  of  field  points  when  it  acts  on  a  sample  of  the 
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odd-shaped  surface's  normal  component  of  velocity.  For 
comparison,  the  corresponding  operator  (3  for  the  case  of  an 
infinite,  continuous  planar  source  can  be  written  as 

a  =  y'M  gh-  9  >  (6.9) 

where  £  is  the  2-D  Fourier  transform  and  G  is  the  Fourier 

N 

transformed  Green's  function  for  planar  Neumann  boundary 

conditions  as  discussed  previously.  The  analogy  between 
the  continuous  and  numerical  processes  can  be  seen  readily: 

The  matrix  [V*]  corresponds  to  the  Fourier  transform,  the 
diagonal  matrix  [Z]  takes  the  place  of  the  multiplication  by 
the  Green's  function,  and  the  matrix  [U]  corresponds  to  the 
Inverse  Fourier  transform.  Notice  that  since  there  are  only 
n  degrees  of  freedom  in  the  odd-shaped  surface  model ,  there 
can  only  be  n  unique  ouput  fields  from  equation  (6.8),  and 
these  fields  are  spanned  by  the  first  n  columns  of  the 
unitary  matrix  [U] .  A  similar  situation  could  arise  in  the 
continuous  case  if  the  boundary  condition  was  known  only  up 
to  a  certain  cutoff  spatial  frequency.  The  propagated  field 
would  then  be  limited  to  the  same  spatial  frequency  range. 

The  reconstruction  of  a  source  plane  field  in  the 
continuous  case  requires  the  application  of  the  inverse  of 
the  operator  <3.  This  operation  consists  of  performing  the 
Fourier  transform  of  the  pressure  in  a  plane  above  the 
source,  dividing  the  transformed  field  by  the  Green's 


function,  and  performing  the  inverse  Fourier  transfrom  on 
the  resulting  field.  If  information  about  the  source  is 
needed  only  up  to  a  finite  spatial  frequency,  then  the 
components  of  higher  spatial  frequencies  produced  by  the 
forward  transform  are  ignored.  Numerical  reconstruction  of 
an  odd-shaped  source  as  given  by  the  generalized  inverse 
of  the  matrix  [GF]  follows  the  same  pattern.  The  pressure 
field  sample  vector  is  decomposed  into  a  superposition  of 
the  columns  of  the  unitary  matrix  [U]  by  multiplication  with 
the  conjugate  transpose  of  [U] .  The  first  n  components  are 
divided  by  the  corresponding  singular  values  of  [£]  while 
the  remaining  m-n  components  are  discarded,  and  the 
resulting  superposition  of  the  columns  of  [V]  is  calculated. 


B . 3  Calculating  the  solution  in  the  presence  of  noise 

The  preceding  discusion  considers  reconstruction  when 
all  quantities  are  known  exactly  in  contrast  with  the  noisy 
data  present  in  experiment.  This  noise  results  from  the 
limitations  on  experimental  measurement  and  from  the 
discrete  approximation  of  continuous  operations  used  to 
process  the  data.  Planar  source  fields  are  reconstructed 
using  the  DFT  of  a  finite  sample  of  the  field  above  the 
source  as  explained  in  chapter  IV.  In  this  case  it  was  shown 
that  the  as  the  spatial  frequency  approximated  by  the  DFT 


grows  higher,  the  noise  level  grows  higher,  and  the 
reciprocal  of  the  Green's  function  needed  in  the 
reconstruction  process  also  grows.  Both  these  factors 
enhance  the  noise  in  the  reconstruction.  The  same  type  of 
problem  can  be  expected  in  any  reconstruction  problem,  and 
methods  for  avoiding  the  problems  with  noise  and  an 
understanding  of  the  effects  of  any  reduction  method  on  the 
solution  are  necessities. 

The  general  approach  to  producing  a  reasonable  solution 
in  the  presence  of  noise  is  to  eliminate  or  reduce  the 
contribution  to  the  solution  from  the  calculated  components 
of  higher  spatial  frequency  modes.  The  relative  magnitude  of 
the  Green's  function  or  singular  values  provides  a  guide  in 
the  noise  reduction  process.  The  use  of  a  filter  function 
was  discussed  in  chapter  VI  for  planar  reconstructions  and  a 
similar  approach  is  suggested  in  the  Geophysics  literature 
for  inversion  problems  using  the  Singular  Value  Decomposi¬ 
tion,  but  no  clear  guidelines  are  provided.48  A  simple 
cutoff  of  any  components  associated  with  a  singular  value 
below  a  chosen  minimum  value  is  another  approach,  and  a 
scheme  for  determining  where  to  set  the  cutoff  based  on  the 
relative  size  of  the  singular  values  and  the  calculated 
component  size  is  presented  in  the  Numerical  Mathematics 

literature.49  The  users  guide43  to  the  Singular  Value 
Decomposition  program  used  in  this  work  maintains  that,  just 


how  to  handle  these  small  singular  values,  is  still  a 
problem  open  to  investigation,  and  future  research  in  this 
area  would  be  valuable  to  a  number  of  disciplines.  Due  to 
the  relatively  small  size  of  the  problem  studied  here  (i.e. 
sixty  degrees  of  freedom)  a  much  simpler  approach  is 
practical . 

The  approach  used  to  generate  the  reconstructions  in 
this  work  is  to  watch  the  reconstructions  develop  as  higher 
and  higher  order  reconstructed  surfaces  modes  are  super¬ 
imposed.  The  process  is  halted  and  reversed  one  step  when 
the  solution  takes  on  a  noisy  appearance.  Research  into 
reducing  the  subjectivity  of  this  step  would  be  desireable. 
However,  notice  that  the  same  cutoff  (or  no  cutoff  at  all) 
is  used  in  the  reconstructions  presented  in  the  final 
section  of  this  chapter.  This  suggests  some  level  of 
consistency  from  experiment  to  experiment.  Furthermore, 
notice  that  the  lower  order  modes  are  not  effected  by  this 
method,  and  these  are  the  modes  which  radiate  most 
efficiently.  In  many  applications,  these  lower  order  modes 
will  be  of  primary  interest,  and  the  neglect  of  tenuous 
information  about  the  highest  order  modes  will  be 
inconsequential . 


The  techniques  of  this  chapter  were  tested  in  a  number 
of  simulated  odd-shaped  source  reconstructions.  A  spherical 
source  boundary  is  represented  in  these  simulations  by  the 
sixty-sided  polyhedron  described  in  the  last  chapter.  A 
simulation  experiment  consists  of  generating  the  theoretical 
pressure  at  a  number  of  points  above  the  odd-shaped  surface, 
thus  forming  a  theoretical  hologram,  and  then  reconstructing 
the  surface  field  from  this  hologram  using  the  techniques 
described  above.  In  this  section  an  appropriate  boundary 
condition  and  a  method  of  calculating  the  field  produced  by 
this  boundary  are  presented. 

C.l.  A  small  piston  set  in  a  rigid  sphere 


A  small  piston  set  in  an  otherwise  rigid  sphere  was 
chosen  as  the  test  case  boundary  condition,  because  it 
provides  a  good  test  of  the  resolving  capabilities  of  the 
reconstruction  technique,  and  the  field  produced  by  such  a 
boundary  is  relatively  easily  calculated.  The  actual 
boundary  condition  used  for  a  sphere  with  radius  a  expressed 
in  spherical  coordinates  is 

f  u  when  0  <  e  <  e 
vn(a,e,0)=  0°  otherwlse 


(6.10)  . 


p- 


The  coordinate  system  is  chosen  to  eliminate  any  dependence 
on  the  angle  0 . 

The  angle  0Q  is  chosen  so  that  the  active  area  of  the 

boundary  is  equal  to  the  area  of  one  triangular  element  in 
the  surface  approximation  model.  This  area  is  one-sixtieth 
of  the  total  sphere's  surface.  Simple  calculation  shows  that 
angle  ©o  satisfies  the  following  equation: 

cos (  e0  )  =  1  -  1/30  (6.11). 

For  this  small  angle,  the  enclosed  surface  deviates  very 
little  from  a  flat  surface,  and  it  is  thus  reasonable  to 
refer  to  this  boundary  section  as  a  small  piston. 

C , 2  Calculating  the  field 


The  pressure  at  any  point  outside  the  rigid  sphere 
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with  a  small  piston  is  given  by  the  infinite  summation: 


P(rf )=  p(rf ,e, 0)= 


ipcu 


h_ ( kr ] 
m 


2  0  SJPm-Jco3*  V]-pm+l[c°s(e0)]}  Pm[cos(e)]  h,(ka) 


m=0 


(6.12) 

where  the  Pm's  are  the  Legendre  polynomials,  the  hm's  are 

the  spherical  Hankel  functions,  and  the  prime  indicates  the 
derivative  with  respect  to  argument.  The  recurrence 
relations  for  these  polynomials  make  them  fast  and  easy  to 


tabulate  numerically.  The  series  also  converges  rapidly  and 
is  approximated  in  this  work  by  the  first  ten  terms. 

D.  Results  from  the  test  case 

In  this  section  the  results  of  simulated  reconstruc¬ 
tions  of  the  piston  set  in  a  rigid  sphere  are  presented.  The 
wavenumber  in  these  simulations  is  4.85  corresponding  to  a 
characteristic  frequency  of  the  sixty-sided  polyhedron 
model,  and  therefore,  the  techniques  of  chapter  V  are 
required  to  form  a  surface  solution  matrix  from  the  over¬ 
determined  system  of  equations  (5.3)  and  (5.5).  A  set  of 
sixty  interior  points  radially  aligned  with  the  surface 
element  centers  but  half  the  distance  from  the  sphere's 
center  were  used  to  generate  an  overdetermined  system  of 
equations  using  the  Interior  Helmholtz  Equation. 

Two  different  hologram  geometries  were  simulated: 
measurements  taken  on  a  planar  grid  of  25x25  microphones 
spaced  symmetrically  0.4X  apart  in  a  plane  0.23X  above  the 
sphere;  and  measurements  taken  at  60  microphones  each 
located  0.23X  radially  above  the  center  of  a  surface  model 
element.  Note  that,  due  to  the  numerical  integration  scheme 
for  determining  the  field  matrix  coefficients,  it  is  best 
not  to  bring  the  field  points  too  close  to  the  surface 
elements.  A  separation  from  any  element  to  any  microphone 


of  half  the  smallest  distance  between  element  centers  should 


be  adequate;  the  distance  between  surface  element  centers  is 
0.48A  in  the  simulations  presented  below. 

D.l  Displaying  the  reconstructed  odd-shaped  surface  fields 

Holographic  reconstructions  generate  large  quantities 
of  information  which  must  be  displayed  graphically  to  be 
useful.  Planar  reconstruction  field  quantities  are  readily 
displayed  by  conventional  3-D  hidden-line  plots  of  the  field 
variable  over  the  2-D  plane,  whereas  odd-shaped  source 
reconstructions  present  a  more  general  graphic  display 
problem.  Novel  graphics  routines  were  developed  to  handle 
this  problem. 

The  technique  for  diplaying  the  reconstructed 
odd-shaped  surface  fields  is  carried  out  in  two  steps. 

First,  enough  imaginary  cuts  are  made  in  the  triangle-faced 
polyhedron  model  of  the  odd-shaped  surface  so  that  it  can  be 
flattened  out  into  a  plane  figure.  The  shape  of  each  surface 
element  triangle  can  be  preserved  in  this  process  provided 
the  original  surface  is  convex.  With  this  process,  some 
elements  of  a  surface  which  is  not  strictly  convex  would 
have  to  be  distorted,  and,  to  avoid  this  distortion,  it 
might  be  preferable  to  process  and  display  sections  of  such 
surfaces  separately.  A  flattened  projection  of  the  sixty 


sided  polyhedron  used  in  this  study  is  presented  in  figure 
6.1.  The  element  numbering  scheme  is  the  same  as  used  in 
figure  5.1  for  the  unprojected  surface. 


Once  a  projection  of  this  type  has  been  achieved,  a 
field  variable  can  be  displayed  by  constructing  a  five-faced 
prism  over  each  triangle  where  the  base  and  top  of  the  prism 
are  identical  triangles  and  the  height  of  each  prism  is 
determined  by  the  field  variable.  A  2-D  projection  of  the 
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resulting  set  of  prisms  can  then  be  generated.  This 
projection  will  display  the  field  variables  over  all  the 
elements  in  a  single  figure  which  conveys  some  idea  about 
their  spatial  distribution  as  well.  The  remaining  figures  of 
this  chapter  serve  as  examples  of  this  projection  scheme. 

D . 2  Reconstruction  from  a  planar  hologram  surface 

Reconstructing  an  odd-shaped  source  from  field  pressure 
measurements  over  a  planar  array  of  microphones  is  of 
special  interest  as  the  Nearfield  Acoustic  Holography  Group 
at  the  Pennsylvania  State  University  possesses  such  a 
measurement  system.  More  generally,  reconstructions  of  this 
type  would  be  experimentally  advantageous  because  a  fixed 
microphone  array  allows  for  very  high  speed  data  acquisition 

as  compared  with  a  scanning  microphone  system. 


The  normal  component  of  intensity  is  an  ideal  quantity 
to  display,  as  it  incorporates  the  normal  component  of 
velocity  and  the  pressure,  and  it  requires  the  proper 
determination  of  phase  between  the  pressure  and  velocity. 
Note  that,  for  comparison  purposes,  the  normal  component  of 
intensity  of  a  small  piston  in  a  rigid  sphere  is  zero 
everywhere  off  the  piston  and  nearly  constant  over  the 
piston . 

The  sixty  sided  polyhedron  is  too  coarse  of  a  model  to 
supply  detailed  quanitative  information  about  the  small 
piston  source,  but  the  following  figures  demonstrate  that  it 
can  be  used  to  determine  the  location  of  the  source  quite 
well,  and  the  intensity  values  produced  are  all  within  3dB 
of  the  theoretical  value  at  the  piston's  center. 

It  is  to  be  expected  that  reconstructions  of  small 
features  will  be  more  difficult  the  further  these  features 
are  from  the  measuring  system.  To  address  this  expectation, 
a  series  of  normal  component  of  intensity  reconstructions 
is  presented  in  which  the  theoretical  piston  is  located 
progressively  further  from  the  hologram  plane. 

Figure  6.2  is  a  3-D  projection  of  the  normal 
component  of  intensity  reconstruction  from  the  theoretical 
pressure  field  of  a  single  piston  centered  on  surface 
element  number  1  [  see  figure  6.1  or  5.1  ].  The  location 
of  the  piston  source  is  clearly  indicated.  A  small  false 
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source  appears  at  element  59  which  is  on  the  opposite  side 
of  the  sphere  from  the  true  source,  however  this  source 
falls  below  7%  the  strength  of  the  true  source.  The  actual 
intensity  value  reconstructed  at  element  1  is  11%  below  the 
theoretical  central  value.  For  this  reconstruction,  out  of 
the  sixty  surface  modes  representable  by  the  surface  model , 
the  modes  associated  with  the  fifty  largest  singular  values 
in  the  decomposition  of  the  field  matrix,  equation  (6.8), 
were  used.  This  holds  true  for  the  next  three  reconstruc¬ 
tions  presented  as  well. 

Figure  6.3  displays  the  reconstructed  normal 
component  of  intensity  for  a  piston  source  centered  on 
surface  element  20.  A  false  source  appears  at  element  53 
which  is  adjacent  to  the  diametrically  opposite  element  of 
element  20.  This  source  is  11%  the  strength  of  the  maximum 
reconstructed  source  strength.  The  reconstructed  source  at 
element  20  is  20%  below  the  theoretical  central  value. 

Figure  6.4  displays  the  reconstructed  normal  component 
of  intensity  for  a  piston  centered  on  on  surface  element  56 
which  is  one  of  the  furthest  elements  from  the  hologram 
plane,  and  it  faces  away  from  the  array.  There  are  several 
small  false  sources  in  the  10%  strength  range,  and  the 
reconstruction  of  the  true  source  is  46%  the  central  theo¬ 


retical  value.  However,  the  existence  of  a  small  localized 
source  at  this  position  is  still  quite  clear. 


6.4  The  reconstructed  normal  component  of  intensity 
a  piston  source  centered  on  surface  element  56. 
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To  see  what  happens  when  an  actual  source  is  poorly 
fit  by  the  surface  model,  a  reconstruction  of  a  piston 
source  which  faces  directly  away  from  the  hologram  plane 
is  presented  in  figure  6.5.  The  orientation  of  the 
theoretical  piston  is  such  that  it  is  centered  on  a 
vertex  of  the  five  bottom-most  elements  and  not  on  any 
single  element.  The  reconstructed  source  strength  at  each  of 
these  five  elements  is  approximately  20%  that  of  the 
theoretical  value  at  the  source's  center.  Again  the 
location  of  the  source  is  clear,  and  it  is  reasonable  to 
claim  that  only  with  a  more  complex  surface  approximation 
could  any  more  detail  about  the  source  be  obtained. 

The  sample  reconstructions  presented  demonstrate  the 
resolution  possible  in  reconstructing  a  high  aspect  ratio 
object,  the  sphere,  from  measurements  over  a  planar  array. 

An  actual  physical  experiment  using  these  techniques,  and 

with  more  sophisticated  surface  approximation  elements,52  is 
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planned  in  the  work  of  another  graduate  student. 

D . 3  Reconstruction  from  a  conforming  hologram  surface 

The  reconstructions  from  a  planar  hologram  presented 
above  required  a  large  number  of  field  measurements  to 
achieve  their  resolution.  Similar  resolution  can  be  achieved 
with  far  fewer  field  measurement  points  if  all  the  points 


fall  on  a  surface  which  is  nearby  and  surrounds  the 
odd-shaped  surface.  As  an  example  of  this,  the  surface  field 
of  the  small  piston  in  a  sphere  is  presented  as 
reconstructed  from  60  field  measurements  taken  over  a 
concentric  sphere  0.23X  away. 

Figure  6.6  displays  the  reconstructed  amplitude  of 
the  normal  component  of  surface  velocity  for  a  piston  source 
centered  on  element  1.  All  sixty  modes  of  the  decomposition 
were  used  in  this  reconstruction  and  those  presented  below; 
these  solutions  can  be  considered  free  of  any  subjective 
filtering.  The  reconstructed  maximum  source  velocity  is  80% 
of  the  theoretical  value.  The  noise  in  this  display  could 
be  decreased  by  filtering  out  those  modes  corresponding  to 
the  smallest  singular  value,  but  it  was  retained  to 
illustrate  that  this  noise  is  not  carried  through  the 
pressure  field  reconstruction  process. 

The  amplitude  of  the  reconstructed  surface  pressure 
field  is  displayed  in  figure  6.7a.  This  surface  field  was 
calculated  from  the  velocity  field  displayed  in  figure  6.6 
and  the  surface  solution  matrix.  The  amplitude  of  the 
theoretical  field  is  shown  in  figure  6.7b  below  the 
experimental  field.  Both  plots  have  the  same  scale  and  the 
agreement  between  the  two  is  quite  clear.  The  apparently 
good  agreement  between  these  two  plots  suggests  that  the 
surface  velocity  field  of  figure  6.6  differs  from  the 


theoretical  field  mostly  by  a  contribution 
from  noisy  but  non-radiating  modes. 

The  reconstructed  normal  component  of  itensity  for  this 
source  is  displayed  in  figure  6.8.  This  plot  clearly  shows 
the  existence  of  a  small  localized  source  and  so  illustrates 
the  accurate  reconstruction  of  phase  information  as  well  as 
amplitude  information. 
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Appendix  A. 

Discrete  Convolution  with  the  DFT 


Discrete  convolutions  such  as  equation  (3.6)  of  chapter 
III,  can  be  performed  exactly  and  quickly  with  the  aid  of 
the  FFT  to  perform  DFT  operations. 

To  introduce  the  DFT  and  IDFT  into  equation  (3.6), 
equation  (3.6)  is  rewritten  with  the  aid  of  the  Kronecker 
Delta  function, 


mn 


-f  1: 

l  0; 


m  =  n 
m  *  n 


The  result  is 


N-l  N—l  N-l  N-l 

tlp.q.z)8  III  I  S'  ( / , m) 3  (r,s,z)5  ,5 

1  =  0  m=0  r=-N  s=-N  D,N  D'  r  ,  p  i  s  ,  q  m 

(  A  .  1  )  . 


The  Delta  function  can  be  written  in  terms  of  the  Fourier 
coefficients  as 


=  i  MS1  e  M(m-n)  where  m=2n. 


mn 


M  ;i=0 

Substituting  this  expansion  for  the  Delta  functions  in 
equation  (A.l)  and  rearranging  the  summations  leaves 


t(p,q,z)= 

M-l  M-l 


i  m-i  n-l  ,  n ,  "N-l  N-l  !n(  ,  , _  , 

—  2  I  e  iflUP+vq)  2  (  1  ,m)  e 

M2  n=0  v—0  [/=0  m=0  D,N 

N- 1  N- 1  _■  rr ,  _  .  _  \ 

2  2  3  ( r ,  s  ,  z )  e  1N{r>1+Sv) 

r=-N  s=-N  0,N 


X 


(A. 2)  . 


Using  the  functions  H»'  and  3'  as  defined  in  equations 


(3.7)  and  (3.12),  the  bracketed  double  sums  in  (A2)  can  be 
rewritten  as  DFT ' s .  By  making  use  of  the  identity 
exp[inr>i/N]  =  exp[  in(  r-2N)  >i/N]  to  arrange  the  second 
bracketed  sum,  equation  (A2)  becomes 


t (p,q,z)= 

M-l  M-l 


1  M-l  M-l  i  tT  ,  _  .  _  . 

—  I  2  e  ^UP+^q) 
M2  >i=0  v=0 


M  1  M-  1  j  tt  /  «  .  „  » 

2  2  l»+nv) 

1  =  0  m=0  D,N 


M-l  M-l 
2  2  G* 

r=0  s=0 


D  ,  N 


(  r ,  s  ,  z )  e  TJ 


-ifllrji+sv) 


(A.  3)  . 


The  summations  in  equation  (A3)  can  now  be  replaced  by  the 
DFT  and  IDFT  operations  as  defined  in  equations  (3.8)  and 
(3.9).  The  final  form  with  these  substitutions  is  given  by 
equations  (3.10)  and  (3.11). 


Appendix  B. 


Quadratic  expansion  and  Integration  of  the  Green's  Functions 

Real-space  Green's  functions 

Rectangular  Elements 

Integrations,  over  rectangular  patches,  of  the  real 
space  Green's  function  are  approximated  by  a  truncated 
Taylor  series  expansion  of  the  function  about  each  patch's 
center  (xQ ,yQ ) : 

G(xq+  x,yQ+y,z-z0)  a:  f(x,y),  with  z-zq  constant,  where 
f(x,y)  =  a  +a  x+o:  v+axy+a  x2+oc  y2+cc  x2y+cx  xy2+a  x2y2 

U  1  2  3  4  3  o  7  8 

and  where  the  are  given  below. 

Integration  of  this  form  over  the  patch  yields 
I«ax  „Ay 

2  2  f ( x , y ) dxdy 

-ax  " -A 
2  2 

=  oc0AxAy  +at  (  Ax)  3  Ay  +a3Ax(  Ay)  3  +gQ ( Ax ) 2 ( Ay ) 2  (B.l). 
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Along  the  z-axis  xQ  =  yQ=  0,  a  different  approach  is 

used:  The  exact  integral  over  a  circular  region  of  radius 

r  around  the  axis  is  performed,  and  the  integration  over  the 
remaining  region,  to  complete  a  square,  is  approximated  by 


expanding  about  the  square's  corner.  In  terms  of  the 
expansion  coefficients,  a^,  a  single  corner  integration 

yields 

PrPr 

f ( x ' , y ' )  dx ' dy 1 = 

0"  ( r2-x ' 2 ) 1 /2 

ao^0r2+(ai+a2^ir3+a3^3r4+(a4+a5)^4r4+(0r6+0t7)^6r5+a^8r 

( B  .  2  ) 

where  x,y,  are  from  the  center  of  the  square,  while  x',y' 
are  from  one  of  its  corners.  The  coefficients  are: 


£o=  2.1460184  x  10_t 
pi=  4.793517  X  10'2 
fi  =  6.2685032  X  10“3 

3 

p4=  1.8252296  X  10"2 

8  =  1.5856291  X  10”3 
6 

pg=  2.889426  x  10"4 


The  exact  integration  over  the 
Neumann  boundary  conditions  is 


2  n 


G  (r',z)  r'dr'de  =  pc' 

N 


central  circular  region 
ikz  _  gik(  r2+z2  ) 1 '/z  j 


for 


( B  .  3  ) 


d  0  d  0 


and  for  Dirichlet  boundary  conditions  it  is 


To  calculate  g^  ^(z)  or  g^  '(z)  from  equation  (3.11), 

S'  (l,m,z)  is  approximated  by  a  polynomial  expansion  as 

given  in  equation  (B.l)  except  when  the  point  (x=0,y=0)  is 
included  in  the  integration  patch  indicated  by  1  and  m. 
Appropriate  values  for  xq  and  y^  are  given  by  x^  and  ym  of 

equation  (3.3)  and  (3.4)  for  g^2^(z),  and  by  xj  =  x^  +  A/2 

and  y}  =  yj  +  A/2  for  gj^(z). 

The  point  (x=0,y=0)  is  centered  in  the  integration 
patch  of  So  N(i=0,m=0,z)  for  g^2^(z).  In  this  case,  the 

approximation  for  G„  (;=0,m=0,z)  is  the  sum  of  the 

appropriate  integral  from  equations  (B.3)  and  (B.4)  with 
r=A/2  plus  four  times  (once  for  each  corner  of  the  patch) 
equation  (B. 2), with  r=A/2  and  with  xq=  yQ=A/2  [to  determine 

the  coefficients  ai(x0'Y0)]- 

When  calculating  g^4)(z)  based  on  equation  (3.11)  and 

0  i  H  v 

the  shifted  sampling  coordinates  x}  and  y^ ,  the  four 

integrations,  which  yield  G„  (7,m,z)  for  (/,m)  =  (0,0), 

( -1 , 0) , ( 0 , -1 ) , ( -1 , -1 ) ,  have  the  point  (x=0,y=0)  at  one 
corner  of  their  integration  patches.  3  (i,m,z)  for  each 

0  /  N 

of  these  four  index  pairs  is  given  by  one-fourth  the 
appropriate  integral  from  equations  (B.3)  and  (B.4)  Plus  the 


corner  integration,  equation  (B.2),  with  r=A  and  with 

x  =y  =A , 

0  1  0 

The  form  gi11^(z)  is  calculated  from  equation  (3.13)  as 

described  except  when  1  -  m  =  0.  For  1  =  m  =  0,  the 
appropriate  integral,  Dirichlet  or  Neumann,  from  equations 

(B.2)  and  (B.4)  is  used  with  r=  (  a2/tt  }  1  ^  . 


The  Expansion  Polynomial  Coefficients  aj.(xl3'Y0)'  where 

a  =  kRQ ,  b  =  kx(j ,  d  =  kyQ ,  k  =  u/c,  and  =  x^  +  y^  +  z: 

are  for  Neumann  Boundary  conditions  (normal  velocity  to 
pressure) : 


“i=  -i-^1b(ia“2-a'3  )ela 


“2=  -i^d(ia-a-a-3)ela 

oc  3  =  -i|^bd(-a_3-3ia“4+3a"s)ela 


Ck4r 


ia  -a 


'+b2(-a"3-3ia“4+3a~5) ]eia 


s=  ~if^k*[ ia~2-a~3+d2 ( -a~3-3ia" 4+3a~ 5 ) ]eia 

.=  -if^-d[-a'3-3ia'4  +  3a's+bz  ( -ia' 4+6a"3  +  15ia'6-15a’ 7  )  ]i 

r=  -if^^b[-a“3-3ia'4  +  3a"5+d2(-ia"4+6a-5  +  15ia'6-15a~7)  ]i 

,=  (-a_3-3ia'4+3a"3) 

+ (b2+d2 ) (-ia-4+6a~5+15ia-6-15a-7) 
+b2d2(a'5+10ia"6-45a'7-105ia"3+105a'3) ]eia  ; 


and  for  Dirichlet  boundary  conditions  (pressure  to 
pressure ) : 


oc0=  £^(a'3-ia'2)ela 

^?b(a"3  +  3ia'4-3a"s  )eia 
a  =  £^d(a'3+3ia"4-3a'3  )eia 

2  c.  TT 

a  =  fe^bd(  ia"4-6a"s-15ia_6  +  15a‘7)ela 

a  +  =  (a-3+3ia"4-3a"3  )+b2  (ia_4-6a_3-15ia_6+15a-7)]eia 

a  =  ^!z[(a-3+3ia-4-3a-3)+d2(ia"4-6a_3-15ia-6+15a'7)]eia 

5  **  TT 

oc  =  ^  [  d  (  ia  4  — 6a  3  — 15ia  ®  +  l5a  7) 

6  nk  TT 

+b2d(-a“3-10ia"6+45a"7+105ia‘s-105a-9 ) ]eia 

a7=  ^~[b(  ia-4-6a~3-15ia~6+15a-7) 

+bd2 (-a"3-10ia“6+45a"7+105ia"3-105a"9 ) ]eia 

a8=  §^[(i  a'4-6a'3-15ia_6  +  15a‘7) 

+2bd(-a_3-10ia-8+45a~7+105ia-8-105a”9) 

+b2d2 (-ia'8+15a_7+105ia"3-420a‘9-945ia'10+945a_1 1 ]e 


Triangular  elements 


Integrations  of  the  Green's  functions  are  done  in  the 
same  manner  as  is  done  for  rectangular  elements,  although 
only  the  first  six  expansion  terms  are  considered.  In  the 


local  coordinate  system  of  each  triangle,  the  vertices  have 
coordinates  (x  ,y  ),  (x2,y2),  and  (X3'Y3)-  The  resuit  of 

integrating  the  Taylor  series  expansion,  in  terms  of  the 

local  vertex  coordinates  and  total  triangle  area  A,54  is 

f (x1 ,y ’ )  dx'dy '= 

A 

a  A+a  A ( x  y  +x  y  +x  v  )+a  A(x2+x2+x2)+oc_A(y2+y2+y2)  . 
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When  calculating  the  diagonal  monopole  coefficient,  equation 
( B  .  3  )  is  evaluated  with  r=(  A/rr )  1 7/3  ,  and  no  expansion  terms 


are  considered. 


Intearated  k-soace  Green's  functions 


The  integrated  k-space  Green's  functions  are  obtained 

by  performing  an  integration  averaging  of  the  analytic 

Green's  function  in  the  area  between  k_  and  k_  where  k_  = 

r  i  r  2  r  i 

(k£  +  k2)l/a  -  Akr  and  kf2=  (k2  +  k2)l/2  +  Akr .  In  this 

study,  Akr=/ZAkx  is  used.  The  resulting  integrated  averages 

for  Neumann  boundary  conditions  are: 

Radiation  region 


(kx,ky,z)  =  - 


2ipckz(elkiz  -eik2z) 

(k2z2  -  k2z2) 
2pck 

” uTTirr" 


k2=k2-k2 

1  r  i 

k2=k2-k2 

2  Yz 


when  z=0 


Mixed  evanescent  and  propagating  region 


9h  ( kx , ky , z )  =  - 


2ipckz(eikiz  -e~k2Z ) 

( k2z2  +  k2z2 ) 

2pck(kA  -ik£) 

(k2  +  k2) 

1  1  2  ' 


k2=k2-k2 

1  r  i 

k2=k2  -k2 

2  r2 


when  z=0 


Evanescent  region 


9h  (kx'ky'z)  =  + 


2ipckz(e~kiz  -e~k2z ) 


( k2z2  -  k2z2) 


-2ipck 
(k.  +  k  ) 


k2=k2  -k2 

1  r  i 

k2=k2  -k2 

2  r  2 


when  z=0 


The  resulting  integrated  averages  for  Dirichlet  boundary 
conditions  are: 


I 
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Appendix  C. 

DFT  Error  in  Approximating  the  Continuous  Fourier  Transform 


The  2-D  Fourier  transform  of  a  function  f(x,y)  is 
given  explicitly  by 


n+°° 


F(kx<k„)= 


o+°° 


f (x,y)e 


-i(kxx+V)dxdy 


(C.l) 


It  is  assumed  that  the  function  f(x,y)  is  reasonably  smooth 
and  is  known  at  a  sufficiently  large  number  of  points 
( xr ' yS )  so  that  one  can  write  f(x,y)=  f(xr,ys)  if  |x-xr|^A/2 

and  |y-yg|£A/2  where  xr=  rA  and  ys=  sa  for  r,s=  0,±l,+2,... 


Equation  (C.l)  can  then  be  written: 


T(kx,k  )=  2  Z  f(x,y) 
/  r  s 


nxr+A/'2  +  2  -i(k  x+k  y) 

e  x  Y  dxdy  ( C . 2 ) . 

xr-A/2  ys~A/2 


The  integral  in  equation  (C.2)  is  readily  evaluated 
yielding : 


4  sin(kxA/2)  sin( kyA/ 2 ) ~ 

k  lc 
x  y 

(  C  .  3  )  . 

The  discrete  Fourier  Transform  of  f(x_,y_)  can  be  written 

r  s 

explicitly  as: 

-i ( k_,x„+k,,y^ ) 

F' (kx,k  )=  a22  Z  f (xr,ys)e  xr  (C.4). 

*  r  s 


F(kx,ky)=Z 


f (xr,yg)e 


-i(kxxr+kyys) 


By  dividing  equation  (C.4)  by  equation  (C.3),  an  estimate  is 
obtained  for  the  error  in  approximating  the  continuous 


Fourier  transform  of  an  approximately  piecewise  constant 
function  f(x,y)  by  the  Discrete  Fourier  transform.  The 


result  is: 


F' (kx,k :  ) 


F ( k„ , k„ ) 


kxA  k  a 


4  sin(k  A/2 )  sin(k  A/2) 
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